first version of multi-slit/pepperpot/what ever aka flexible collimator....

first version of multi-slit/pepperpot/what ever aka flexible collimator. However way to slow to be practical

src/Algorithms/LieMapper.cpp

@@ -25,6 +25,7 @@
 #include "AbsBeamline/Degrader.h"
 #include "AbsBeamline/Diagnostic.h"
 #include "AbsBeamline/Drift.h"
+#include "AbsBeamline/FlexibleCollimator.h"
 #include "AbsBeamline/ElementBase.h"
 #include "AbsBeamline/Lambertson.h"
 #include "AbsBeamline/Marker.h"
@@ -157,6 +158,9 @@ void LieMapper::visitDrift(const Drift &drift) {
     applyDrift(flip_s * drift.getElementLength());
 }
 
+void LieMapper::visitFlexibleCollimator(const FlexibleCollimator &coll) {
+    applyDrift(flip_s * coll.getElementLength());
+}
 
 void LieMapper::visitLambertson(const Lambertson &lamb) {
     // Assume the particle go through the magnet's window.

src/Algorithms/LieMapper.h

@@ -111,6 +111,9 @@ public:
     /// Apply the algorithm to a Drift.
     virtual void visitDrift(const Drift &);
 
+    /// Apply the algorithm to a flexible collimator
+    virtual void visitFlexibleCollimator(const FlexibleCollimator &);
+
     /// Apply the algorithm to a Lambertson.
     virtual void visitLambertson(const Lambertson &);
 

src/Algorithms/NilTracker.h

@@ -25,6 +25,7 @@ class Degrader;
 class Diagnostic;
 class Drift;
 class ElementBase;
+class FlexibleCollimator;
 class Lambertson;
 class Marker;
 class Monitor;
@@ -60,6 +61,7 @@ public:
     NIL_VISITELEMENT(Degrader)
     NIL_VISITELEMENT(Diagnostic)
     NIL_VISITELEMENT(Drift)
+    NIL_VISITELEMENT(FlexibleCollimator)
     NIL_VISITELEMENT(Lambertson)
     NIL_VISITELEMENT(Marker)
     NIL_VISITELEMENT(Monitor)

src/Algorithms/ParallelCyclotronTracker.cpp

@@ -693,6 +693,15 @@ void ParallelCyclotronTracker::visitDrift(const Drift &drift) {
     myElements.push_back(dynamic_cast<Drift *>(drift.clone()));
 }
 
+/**
+ *  
+ *
+ *  @param 
+ */
+void ParallelCyclotronTracker::visitFlexibleCollimator(const FlexibleCollimator &) {
+
+}
+
 /**
  *
  *
@@ -3573,4 +3582,4 @@ void ParallelCyclotronTracker::injectBunch_m(bool& flagTransition) {
         // After this, numBunch_m is wrong but not needed anymore...
         numBunch_m--;
     }
-}
\ No newline at end of file
+}

src/Algorithms/ParallelCyclotronTracker.h

@@ -117,6 +117,9 @@ public:
     /// Apply the algorithm to a Drift.
     virtual void visitDrift(const Drift &);
 
+    /// Apply the algorithm to a flexible collimator
+    virtual void visitFlexibleCollimator(const FlexibleCollimator &);
+
     /// Apply the algorithm to a Lambertson.
     virtual void visitLambertson(const Lambertson &);
 

src/Algorithms/ParallelSliceTracker.h

@@ -67,6 +67,7 @@ public:
     virtual void visitDegrader(const Degrader &);
     virtual void visitDiagnostic(const Diagnostic &);
     virtual void visitDrift(const Drift &);
+    virtual void visitFlexibleCollimator(const FlexibleCollimator &);
     virtual void visitLambertson(const Lambertson &);
     virtual void visitMarker(const Marker &);
     virtual void visitMonitor(const Monitor &);
@@ -337,4 +338,4 @@ inline void ParallelSliceTracker::writeLastStepPhaseSpace(const long long step,
         INFOMSG("* Invalid bunch! No statistics dumped." << endl);
 }
 
-#endif
\ No newline at end of file
+#endif

src/Algorithms/ParallelTTracker.cpp

@@ -722,6 +722,7 @@ void ParallelTTracker::computeWakefield(IndexMap::value_t &elements) {
 
             for (unsigned int i = 0; i < localNum; ++ i) {
                 itsBunch_m->R[i] = referenceToBeamCSTrafo.transformTo(itsBunch_m->R[i]);
+                itsBunch_m->P[i] = referenceToBeamCSTrafo.rotateTo(itsBunch_m->P[i]);
                 itsBunch_m->Ef[i] = referenceToBeamCSTrafo.rotateTo(itsBunch_m->Ef[i]);
             }
 
@@ -729,6 +730,7 @@ void ParallelTTracker::computeWakefield(IndexMap::value_t &elements) {
 
             for (unsigned int i = 0; i < localNum; ++ i) {
                 itsBunch_m->R[i] = beamToReferenceCSTrafo.transformTo(itsBunch_m->R[i]);
+                itsBunch_m->P[i] = beamToReferenceCSTrafo.rotateTo(itsBunch_m->P[i]);
                 itsBunch_m->Ef[i] = beamToReferenceCSTrafo.rotateTo(itsBunch_m->Ef[i]);
             }
 

src/Algorithms/ParallelTTracker.h

@@ -36,6 +36,7 @@
 #include "AbsBeamline/Diagnostic.h"
 #include "AbsBeamline/Degrader.h"
 #include "AbsBeamline/Drift.h"
+#include "AbsBeamline/FlexibleCollimator.h"
 #include "AbsBeamline/ElementBase.h"
 #include "AbsBeamline/Lambertson.h"
 #include "AbsBeamline/Marker.h"
@@ -120,6 +121,9 @@ public:
     /// Apply the algorithm to a Drift.
     virtual void visitDrift(const Drift &);
 
+    /// Apply the algorithm to a flexible collimator
+    virtual void visitFlexibleCollimator(const FlexibleCollimator &);
+
     /// Apply the algorithm to a Lambertson.
     virtual void visitLambertson(const Lambertson &);
 
@@ -353,6 +357,11 @@ inline void ParallelTTracker::visitDrift(const Drift &drift) {
 }
 
 
+inline void ParallelTTracker::visitFlexibleCollimator(const FlexibleCollimator &coll) {
+    itsOpalBeamline_m.visit(coll, *this, itsBunch_m);
+}
+
+
 inline void ParallelTTracker::visitLambertson(const Lambertson &lamb) {
     itsOpalBeamline_m.visit(lamb, *this, itsBunch_m);
 }

src/Algorithms/StatisticalErrors.h

@@ -31,6 +31,7 @@ class Degrader;
 class Diagnostic;
 class Drift;
 class ElementBase;
+class FlexibleCollimator;
 class Lambertson;
 class Marker;
 class Monitor;
@@ -72,6 +73,7 @@ public:
     SE_VISITELEMENT(Degrader)
     SE_VISITELEMENT(Diagnostic)
     SE_VISITELEMENT(Drift)
+    SE_VISITELEMENT(FlexibleCollimator)
     SE_VISITELEMENT(Lambertson)
     SE_VISITELEMENT(Marker)
     SE_VISITELEMENT(Monitor)

src/Algorithms/ThickMapper.cpp

@@ -26,6 +26,7 @@
 #include "AbsBeamline/Diagnostic.h"
 #include "AbsBeamline/Drift.h"
 #include "AbsBeamline/ElementBase.h"
+#include "AbsBeamline/FlexibleCollimator.h"
 #include "AbsBeamline/Lambertson.h"
 #include "AbsBeamline/Marker.h"
 #include "AbsBeamline/Monitor.h"
@@ -128,6 +129,9 @@ void ThickMapper::visitDrift(const Drift &drift) {
     applyDrift(flip_s * drift.getElementLength());
 }
 
+void ThickMapper::visitFlexibleCollimator(const FlexibleCollimator &coll) {
+    applyDrift(flip_s * coll.getElementLength());
+}
 
 void ThickMapper::visitLambertson(const Lambertson &lamb) {
     // Assume the particle go through the magnet's window.

src/Algorithms/ThickMapper.h

@@ -113,6 +113,9 @@ public:
     /// Apply the algorithm to a Drift.
     virtual void visitDrift(const Drift &);
 
+    /// Apply the algorithm to a flexible collimator
+    virtual void visitFlexibleCollimator(const FlexibleCollimator &);
+
     /// Apply the algorithm to a Lambertson.
     virtual void visitLambertson(const Lambertson &);
 

src/Algorithms/ThickTracker.h

@@ -20,7 +20,7 @@
 
 #include "Algorithms/Tracker.h"
 #include "Steppers/BorisPusher.h"
-#include "Structure/DataSink.h"                                                                      
+#include "Structure/DataSink.h"
 
 #include "Algorithms/OrbitThreader.h"
 #include "Algorithms/IndexMap.h"
@@ -31,6 +31,7 @@
 #include "AbsBeamline/Diagnostic.h"
 #include "AbsBeamline/Degrader.h"
 #include "AbsBeamline/Drift.h"
+#include "AbsBeamline/FlexibleCollimator.h"
 #include "AbsBeamline/ElementBase.h"
 #include "AbsBeamline/Lambertson.h"
 #include "AbsBeamline/Marker.h"
@@ -48,7 +49,7 @@
 #include "AbsBeamline/Solenoid.h"
 #include "AbsBeamline/ParallelPlate.h"
 #include "AbsBeamline/CyclotronValley.h"
-                       
+
 #include "Elements/OpalBeamline.h"
 
 class BMultipoleField;
@@ -118,13 +119,13 @@ public:
     //  If [b]revTrack[/b] is true, we track against the beam.
     explicit ThickTracker(const Beamline &bl, PartBunchBase<double, 3> *bunch,
 			  DataSink &ds,
-                          const PartData &data, 
+                          const PartData &data,
 			  bool revBeam, bool revTrack,
 			  const std::vector<unsigned long long> &maxSTEPS,
 			  double zstart,
 			  const std::vector<double> &zstop,
 			  const std::vector<double> &dt);
-    
+
     virtual ~ThickTracker();
 
 
@@ -148,6 +149,9 @@ public:
     /// Apply the algorithm to a Drift.
     virtual void visitDrift(const Drift &);
 
+    /// Apply the algorithm to a flexible collimator
+    virtual void visitFlexibleCollimator(const FlexibleCollimator &);
+
     /// Apply the algorithm to a Lambertson.
     virtual void visitLambertson(const Lambertson &);
 
@@ -204,7 +208,7 @@ public:
 
     void updateReferenceParticle(const BorisPusher &pusher);
     void updateRFElement(std::string elName, double maxPhase);
-    void prepareSections();  
+    void prepareSections();
     void saveCavityPhases();
     void restoreCavityPhases();
     void changeDT();
@@ -225,10 +229,10 @@ private:
     void applyExitFringe(double edge, double curve,
                          const BMultipoleField &field, double scale);
 
-    Vector_t RefPartR_m; 
-    Vector_t RefPartP_m; 
+    Vector_t RefPartR_m;
+    Vector_t RefPartP_m;
 
-    DataSink *itsDataSink_m; 
+    DataSink *itsDataSink_m;
 
     OpalBeamline itsOpalBeamline_m;
 
@@ -241,10 +245,10 @@ private:
     /// where to stop
     std::queue<double> zStop_m;
 
-    double dtCurrentTrack_m;  
+    double dtCurrentTrack_m;
     std::queue<double> dtAllTracks_m;
 
-    /// The maximal number of steps the system is integrated per TRACK 
+    /// The maximal number of steps the system is integrated per TRACK
     std::queue<unsigned long long> localTrackSteps_m;
 
     CoordinateSystemTrafo referenceToLabCSTrafo_m;
@@ -287,6 +291,11 @@ inline void ThickTracker::visitDrift(const Drift &drift) {
 }
 
 
+inline void ThickTracker::visitFlexibleCollimator(const FlexibleCollimator &coll) {
+    itsOpalBeamline_m.visit(coll, *this, itsBunch_m);
+}
+
+
 inline void ThickTracker::visitLambertson(const Lambertson &lamb) {
     itsOpalBeamline_m.visit(lamb, *this, itsBunch_m);
 }
@@ -357,4 +366,4 @@ inline void ThickTracker::visitCyclotronValley(const CyclotronValley &cv) {
     itsOpalBeamline_m.visit(cv, *this, itsBunch_m);
 }
 
-#endif // OPAL_ThickTracker_HH
+#endif // OPAL_ThickTracker_HH
\ No newline at end of file

src/Algorithms/TransportMapper.cpp

@@ -26,6 +26,7 @@
 #include "AbsBeamline/Degrader.h"
 #include "AbsBeamline/Drift.h"
 #include "AbsBeamline/ElementBase.h"
+#include "AbsBeamline/FlexibleCollimator.h"
 #include "AbsBeamline/Lambertson.h"
 #include "AbsBeamline/Monitor.h"
 #include "AbsBeamline/Multipole.h"
@@ -160,6 +161,9 @@ void TransportMapper::visitDrift(const Drift &drift) {
     applyDrift(flip_s * drift.getElementLength());
 }
 
+void TransportMapper::visitFlexibleCollimator(const FlexibleCollimator &coll) {
+    applyDrift(flip_s * coll.getElementLength());
+}
 
 void TransportMapper::visitLambertson(const Lambertson &lamb) {
     // Assume the particle go through the magnet's window.

src/Algorithms/TransportMapper.h

@@ -121,6 +121,9 @@ public:
     /// Apply the algorithm to a Drift.
     virtual void visitDrift(const Drift &);
 
+    /// Apply the algorithm to a flexible collimator
+    virtual void visitFlexibleCollimator(const FlexibleCollimator &);
+
     /// Apply the algorithm to a Lambertson.
     virtual void visitLambertson(const Lambertson &);
 

src/Classic/AbsBeamline/BeamlineVisitor.h

@@ -41,6 +41,7 @@ class Corrector;
 class Degrader;
 class Diagnostic;
 class Drift;
+class FlexibleCollimator;
 class Lambertson;
 class Marker;
 class Monitor;
@@ -125,6 +126,9 @@ public:
     /// Apply the algorithm to a drift space.
     virtual void visitDrift(const Drift &) = 0;
 
+    /// Apply the algorithm to a flexible collimator
+    virtual void visitFlexibleCollimator(const FlexibleCollimator &) = 0;
+
     /// Apply the algorithm to an Ring
     virtual void visitRing(const Ring &) = 0;
 

src/Classic/AbsBeamline/CMakeLists.txt

@@ -17,6 +17,7 @@ set (_SRCS
   Drift.cpp
   ElementBase.cpp
   ElementImage.cpp
+  FlexibleCollimator.cpp
   Integrator.cpp
   Lambertson.cpp
   Marker.cpp
@@ -67,6 +68,7 @@ set (HDRS
     Drift.h
     ElementBase.h
     ElementImage.h
+    FlexibleCollimator.h
     Integrator.h
     Lambertson.h
     Marker.h
@@ -96,4 +98,4 @@ set (HDRS
     VariableRFCavity.h
 )
 
-install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline")
+install (FILES ${HDRS} DESTINATION "${CMAKE_INSTALL_PREFIX}/include/AbsBeamline")
\ No newline at end of file

src/Classic/AbsBeamline/FlexibleCollimator.cpp

+#include "AbsBeamline/FlexibleCollimator.h"
+#include "Physics/Physics.h"
+#include "Algorithms/PartBunchBase.h"
+#include "AbsBeamline/BeamlineVisitor.h"
+#include "Fields/Fieldmap.h"
+#include "Structure/LossDataSink.h"
+#include "Utilities/Options.h"
+#include "Solvers/ParticleMatterInteractionHandler.hh"
+#include "Utilities/Util.h"
+
+#include <memory>
+
+extern Inform *gmsg;
+
+// Class FlexibleCollimator
+// ------------------------------------------------------------------------
+
+FlexibleCollimator::FlexibleCollimator():
+    Component(),
+    filename_m(""),
+    informed_m(false),
+    losses_m(0),
+    losses1_m(0),
+    losses2_m(0),
+    minR_m(0),
+    maxR_m(0),
+    lossDs_m(nullptr),
+    parmatint_m(NULL){}
+
+
+FlexibleCollimator::FlexibleCollimator(const FlexibleCollimator &right):
+    Component(right),
+    holes_m(right.holes_m.begin(), right.holes_m.end()),
+    bb_m(right.bb_m),
+    filename_m(right.filename_m),
+    informed_m(right.informed_m),
+    losses_m(0),
+    losses1_m(0),
+    losses2_m(0),
+    minR_m(0),
+    maxR_m(0),
+    lossDs_m(nullptr),
+    parmatint_m(NULL) {
+}
+
+
+FlexibleCollimator::FlexibleCollimator(const std::string &name):
+    Component(name),
+    filename_m(""),
+    informed_m(false),
+    losses_m(0),
+    losses1_m(0),
+    losses2_m(0),
+    minR_m(0),
+    maxR_m(0),
+    lossDs_m(nullptr),
+    parmatint_m(NULL)
+{}
+
+
+FlexibleCollimator::~FlexibleCollimator() {
+    if (online_m)
+        goOffline();
+}
+
+
+void FlexibleCollimator::accept(BeamlineVisitor &visitor) const {
+    visitor.visitFlexibleCollimator(*this);
+}
+
+
+bool FlexibleCollimator::isStopped(const Vector_t &R, const Vector_t &P, double recpgamma) {
+    const double z = R(2);
+
+    if ((z < 0.0) ||
+        (z > getElementLength()) ||
+        (!isInsideTransverse(R))) return false;
+
+    if (!bb_m.isInside(R)) {
+        ++ losses1_m;
+        if (R[0] < minR_m[0]) minR_m[0] = R[0];
+        else if (R[0] > maxR_m[0]) maxR_m[0] = R[0];
+        if (R[1] < minR_m[1]) minR_m[1] = R[1];
+        else if (R[1] > maxR_m[1]) maxR_m[1] = R[1];
+        return true;
+    }
+
+    for (mslang::Base *func: holes_m) {
+        if (func->isInside(R)) return false;
+    }
+
+    ++ losses2_m;
+    return true;
+}
+
+bool FlexibleCollimator::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
+    const Vector_t &R = RefPartBunch_m->R[i];
+    const Vector_t &P = RefPartBunch_m->P[i];
+    const double &dt = RefPartBunch_m->dt[i];
+    const double recpgamma = Physics::c * dt / sqrt(1.0 + dot(P, P));
+    bool pdead = isStopped(R, P, recpgamma);
+
+    if (pdead) {
+        if (lossDs_m) {
+            double frac = -R(2) / P(2) * recpgamma;
+            lossDs_m->addParticle(R, P,
+                                  RefPartBunch_m->ID[i],
+                                  t + frac * dt, 0);
+        }
+        ++ losses_m;
+    }
+
+    return pdead;
+}
+
+bool FlexibleCollimator::applyToReferenceParticle(const Vector_t &R, const Vector_t &P, const double &t, Vector_t &E, Vector_t &B) {
+    *gmsg << __DBGMSG__
+          << std::setw(10) << losses1_m
+          << std::setw(10) << losses2_m
+          << std::setw(18) << minR_m[0]
+          << std::setw(18) << minR_m[1]
+          << std::setw(18) << maxR_m[0]
+          << std::setw(18) << maxR_m[1]
+          << std::setw(18) << bb_m.center_m[0]
+          << std::setw(18) << bb_m.center_m[1]
+          << std::setw(18) << bb_m.width_m
+          << std::setw(18) << bb_m.height_m
+          << endl;
+    losses1_m = losses2_m = 0;
+    minR_m = Vector_t(0.0);
+    maxR_m = Vector_t(0.0);
+
+    return false;
+}
+
+// rectangle collimators in cyclotron cyclindral coordiantes
+// without particlematterinteraction, the particle hitting collimator is deleted directly
+bool FlexibleCollimator::checkCollimator(PartBunchBase<double, 3> *bunch, const int turnnumber, const double t, const double tstep) {
+
+    return false;
+}
+
+void FlexibleCollimator::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) {
+    RefPartBunch_m = bunch;
+    endField = startField + getElementLength();
+
+    parmatint_m = getParticleMatterInteraction();
+
+    // if (!parmatint_m) {
+    if (filename_m == std::string(""))
+        lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
+    else
+        lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
+    // }
+
+    goOnline(-1e6);
+}
+
+void FlexibleCollimator::initialise(PartBunchBase<double, 3> *bunch) {
+    RefPartBunch_m = bunch;
+
+    parmatint_m = getParticleMatterInteraction();
+
+    // if (!parmatint_m) {
+    if (filename_m == std::string(""))
+        lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(getName(), !Options::asciidump));
+    else
+        lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m.substr(0, filename_m.rfind(".")), !Options::asciidump));
+    // }
+
+    goOnline(-1e6);
+}
+
+void FlexibleCollimator::finalise()
+{
+    if (online_m)
+        goOffline();
+    *gmsg << "* Finalize probe" << endl;
+}
+
+void FlexibleCollimator::goOnline(const double &) {
+    print();
+
+    online_m = true;
+}
+
+void FlexibleCollimator::print() {
+    if (RefPartBunch_m == NULL) {
+        if (!informed_m) {
+            std::string errormsg = Fieldmap::typeset_msg("BUNCH SIZE NOT SET", "warning");
+            ERRORMSG(errormsg << endl);
+            if (Ippl::myNode() == 0) {
+                std::ofstream omsg("errormsg.txt", std::ios_base::app);
+                omsg << errormsg << std::endl;
+                omsg.close();
+            }
+            informed_m = true;
+        }
+        return;
+    }
+
+    *gmsg << level3;
+}
+
+void FlexibleCollimator::goOffline() {
+    if (online_m && lossDs_m)
+        lossDs_m->save();
+    lossDs_m.reset(0);
+    online_m = false;
+}
+
+void FlexibleCollimator::getDimensions(double &zBegin, double &zEnd) const {
+    zBegin = 0.0;
+    zEnd = getElementLength();
+}
+
+ElementBase::ElementType FlexibleCollimator::getType() const {
+    return COLLIMATOR;
+}
+
+void FlexibleCollimator::setDescription(const std::string &desc) {
+    mslang::Function *fun;
+
+    if (!mslang::parse(desc, fun))
+        throw GeneralClassicException("FlexibleCollimator::setDescription",
+                                      "Couldn't parse input file");
+
+    fun->apply(holes_m);
+
+    if (holes_m.size() == 0) return;
+
+    for (auto it: holes_m) {
+        it->computeBoundingBox();
+    }
+
+    mslang::Base *first = holes_m.front();
+    const mslang::BoundingBox &bb = first->bb_m;
+
+    Vector_t llc(bb.center_m[0] - 0.5 * bb.width_m,
+                 bb.center_m[1] - 0.5 * bb.height_m,
+                 0.0);
+    Vector_t urc(bb.center_m[0] + 0.5 * bb.width_m,
+                 bb.center_m[1] + 0.5 * bb.height_m,
+                 0.0);
+
+    for (auto it: holes_m) {
+        const mslang::BoundingBox &bb = it->bb_m;
+        llc[0] = std::min(llc[0], bb.center_m[0] - 0.5 * bb.width_m);
+        llc[1] = std::min(llc[1], bb.center_m[1] - 0.5 * bb.height_m);
+        urc[0] = std::max(urc[0], bb.center_m[0] + 0.5 * bb.width_m);
+        urc[1] = std::max(urc[1], bb.center_m[1] + 0.5 * bb.height_m);