CCollimator.cpp 5.83 KB
Newer Older
kraus's avatar
kraus committed
1
#include "AbsBeamline/CCollimator.h"
2

kraus's avatar
kraus committed
3
#include "AbsBeamline/BeamlineVisitor.h"
4
#include "Algorithms/PartBunchBase.h"
kraus's avatar
kraus committed
5
#include "Fields/Fieldmap.h"
6
#include "Solvers/ParticleMatterInteractionHandler.hh"
kraus's avatar
kraus committed
7 8
#include "Structure/LossDataSink.h"

9 10
#include <cmath>
#include <fstream>
kraus's avatar
kraus committed
11 12 13

extern Inform *gmsg;

14 15
CCollimator::CCollimator():CCollimator("")
{}
kraus's avatar
kraus committed
16

17 18
CCollimator::CCollimator(const std::string &name):
    PluginElement(name) {
19
    setDimensions(0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0);
20
    setGeom(0.0);
21
}
kraus's avatar
kraus committed
22 23

CCollimator::CCollimator(const CCollimator &right):
24 25
    PluginElement(right),
    informed_m(right.informed_m) {
26 27 28 29
    setDimensions(right.xstart_m, right.xend_m,
                  right.ystart_m, right.yend_m,
                  right.zstart_m, right.zend_m,
                  right.width_m);
30
    setGeom(width_m);
kraus's avatar
kraus committed
31 32
}

33
CCollimator::~CCollimator() {}
kraus's avatar
kraus committed
34 35 36 37 38

void CCollimator::accept(BeamlineVisitor &visitor) const {
    visitor.visitCCollimator(*this);
}

39
bool CCollimator::doPreCheck(PartBunchBase<double, 3> *bunch) {
kraus's avatar
kraus committed
40 41 42 43
    Vector_t rmin, rmax;
    bunch->get_bounds(rmin, rmax);

    if (rmax(2) >= zstart_m && rmin(2) <= zend_m) {
44 45 46 47 48 49 50 51 52
        // interested in absolute minimum and maximum
        double xmin = std::min(std::abs(rmin(0)), std::abs(rmax(0)));
        double xmax = std::max(std::abs(rmin(0)), std::abs(rmax(0)));
        double ymin = std::min(std::abs(rmin(1)), std::abs(rmax(1)));
        double ymax = std::max(std::abs(rmin(1)), std::abs(rmax(1)));
        double rbunch_min = std::hypot(xmin, ymin);
        double rbunch_max = std::hypot(xmax, ymax);

        if( rbunch_max > rmin_m && rbunch_min < rmax_m ){ // check similar to z
53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76
            return true;
        }
    }
    return false;
}

// rectangle collimators in cyclotron cylindrical coordinates
// when there is no particlematterinteraction, the particle hitting collimator is deleted directly
bool CCollimator::doCheck(PartBunchBase<double, 3> *bunch, const int /*turnnumber*/, const double /*t*/, const double /*tstep*/) {

    bool flagNeedUpdate = false;
    size_t tempnum = bunch->getLocalNum();
    int pflag = 0;
    // now check each particle in bunch
    for (unsigned int i = 0; i < tempnum; ++i) {
        if (bunch->PType[i] == ParticleType::REGULAR && bunch->R[i](2) < zend_m && bunch->R[i](2) > zstart_m ) {
            // only now careful check in r
            pflag = checkPoint(bunch->R[i](0), bunch->R[i](1));
            /// bunch->Bin[i] != -1 makes sure the particle is not stored in more than one collimator
            if ((pflag != 0) && (bunch->Bin[i] != -1)) {
                if (!parmatint_m)
                    lossDs_m->addParticle(bunch->R[i], bunch->P[i], bunch->ID[i]);
                bunch->Bin[i] = -1;
                flagNeedUpdate = true;
kraus's avatar
kraus committed
77 78 79
            }
        }
    }
80 81 82 83
    return flagNeedUpdate;
}

bool CCollimator::doFinaliseCheck(PartBunchBase<double, 3> *bunch, bool flagNeedUpdate) {
kraus's avatar
kraus committed
84 85
    reduce(&flagNeedUpdate, &flagNeedUpdate + 1, &flagNeedUpdate, OpBitwiseOrAssign());
    if (flagNeedUpdate && parmatint_m) {
86 87
        Vector_t rmin, rmax;
        bunch->get_bounds(rmin, rmax);
88 89 90
        std::pair<Vector_t, double> boundingSphere;
        boundingSphere.first = 0.5 * (rmax + rmin);
        boundingSphere.second = euclidean_norm(rmax - boundingSphere.first);
kraus's avatar
kraus committed
91 92 93 94 95
        parmatint_m->apply(bunch, boundingSphere);
    }
    return flagNeedUpdate;
}

96
void CCollimator::doInitialise(PartBunchBase<double, 3> *bunch) {
kraus's avatar
kraus committed
97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112
    parmatint_m = getParticleMatterInteraction();
}

void CCollimator::goOnline(const double &) {
    print();
    // PosX_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    // PosY_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    // PosZ_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    // MomentumX_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    // MomentumY_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    // MomentumZ_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    // time_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    // id_m.reserve((int)(1.1 * RefPartBunch_m->getLocalNum()));
    online_m = true;
}

113 114 115 116
void CCollimator::doFinalise() {
    *gmsg << "* Finalize cyclotron collimator " << getName() << endl;
}

kraus's avatar
kraus committed
117
void CCollimator::print() {
118
    if (RefPartBunch_m == nullptr) {
kraus's avatar
kraus committed
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135
        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 << "* CCollimator angle start " << xstart_m << " (Deg) angle end " << ystart_m << " (Deg) "
          << "R start " << xend_m << " (mm) R rend " << yend_m << " (mm)" << endl;
}

136
void CCollimator::setDimensions(double xstart, double xend, double ystart, double yend, double zstart, double zend, double width) {
137
    setDimensions(xstart, xend, ystart, yend);
138 139 140
    zstart_m = zstart;
    zend_m   = zend;
    width_m  = width;
141 142 143 144
    // zstart and zend are independent from x, y
    if (zstart_m > zend_m){
        std::swap(zstart_m, zend_m);
    }
145
    setGeom(width_m);
kraus's avatar
kraus committed
146 147 148 149 150 151 152 153 154 155 156
}

void CCollimator::getDimensions(double &zBegin, double &zEnd) const {
    zBegin = 0.0;
    zEnd = getElementLength();
}

ElementBase::ElementType CCollimator::getType() const {
    return CCOLLIMATOR;
}

157
void CCollimator::doSetGeom() {
158
    // calculate maximum r from these
159 160 161 162
    rmax_m = -std::numeric_limits<double>::max();
    for (int i=0; i<4; i++) {
        double r = std::hypot(geom_m[i].x, geom_m[i].y);
        rmax_m = std::max(rmax_m, r);
kraus's avatar
kraus committed
163
    }
164 165 166 167 168
    // some debug printout
    // for (int i=0; i<4; i++) {
    //     *gmsg << "point " << i << " ( " << geom_m[i].x << ", " << geom_m[i].y << ")" << endl;
    // }
    // *gmsg << "rmin " << rmin_m << " rmax " << rmax_m << endl;
169
}