Monitor.cpp 7.21 KB
Newer Older
gsell's avatar
gsell committed
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
// ------------------------------------------------------------------------
// $RCSfile: Monitor.cpp,v $
// ------------------------------------------------------------------------
// $Revision: 1.1.1.1 $
// ------------------------------------------------------------------------
// Copyright: see Copyright.readme
// ------------------------------------------------------------------------
//
// Class: Monitor
//   Defines the abstract interface for a beam position monitor.
//
// ------------------------------------------------------------------------
// Class category: AbsBeamline
// ------------------------------------------------------------------------
//
// $Date: 2000/03/27 09:32:31 $
// $Author: fci $
//
// ------------------------------------------------------------------------
kraus's avatar
kraus committed
20
#include "AbsBeamline/Monitor.h"
gsell's avatar
gsell committed
21
#include "Physics/Physics.h"
22
#include "Algorithms/PartBunchBase.h"
gsell's avatar
gsell committed
23
#include "AbsBeamline/BeamlineVisitor.h"
24
#include "Fields/Fieldmap.h"
adelmann's avatar
adelmann committed
25
#include "Structure/LossDataSink.h"
adelmann's avatar
adelmann committed
26
#include "Utilities/Options.h"
27 28 29
#include "Utilities/Util.h"
#include <boost/filesystem.hpp>
#include "AbstractObjects/OpalData.h"
30
#include "Structure/MonitorStatisticsWriter.h"
adelmann's avatar
adelmann committed
31

gsell's avatar
gsell committed
32
#include <fstream>
33
#include <memory>
gsell's avatar
gsell committed
34 35 36 37


// Class Monitor
// ------------------------------------------------------------------------
38 39 40 41

std::map<double, SetStatistics> Monitor::statFileEntries_sm;
const double Monitor::halfLength_s = 0.005;

gsell's avatar
gsell committed
42
Monitor::Monitor():
43
    Monitor("")
gsell's avatar
gsell committed
44 45 46 47 48 49 50
{}


Monitor::Monitor(const Monitor &right):
    Component(right),
    filename_m(right.filename_m),
    plane_m(right.plane_m),
kraus's avatar
kraus committed
51
    type_m(right.type_m),
52
    numPassages_m(0)
gsell's avatar
gsell committed
53 54 55
{}


56
Monitor::Monitor(const std::string &name):
gsell's avatar
gsell committed
57 58 59
    Component(name),
    filename_m(""),
    plane_m(OFF),
kraus's avatar
kraus committed
60
    type_m(SPATIAL),
61
    numPassages_m(0)
gsell's avatar
gsell committed
62 63 64 65 66 67 68 69 70 71 72
{}


Monitor::~Monitor()
{}


void Monitor::accept(BeamlineVisitor &visitor) const {
    visitor.visitMonitor(*this);
}

73
bool Monitor::apply(const size_t &i, const double &t, Vector_t &/*E*/, Vector_t &/*B*/) {
gsell's avatar
gsell committed
74 75
    const Vector_t &R = RefPartBunch_m->R[i];
    const Vector_t &P = RefPartBunch_m->P[i];
76 77 78 79
    const double &dt = RefPartBunch_m->dt[i];
    const double recpgamma = Physics::c * dt / Util::getGamma(P);
    const double middle = 0.5 * getElementLength();
    if (online_m && type_m == SPATIAL) {
80 81
        if (dt * R(2) < dt * middle && // multiply with dt to allow back tracking
            dt * (R(2) + P(2) * recpgamma) > dt * middle) {
82 83 84 85
            double frac = (middle - R(2)) / (P(2) * recpgamma);

            lossDs_m->addParticle(R + frac * recpgamma * P,
                                  P, RefPartBunch_m->ID[i], t + frac * dt, 0);
kraus's avatar
kraus committed
86
        }
gsell's avatar
gsell committed
87 88 89 90 91
    }

    return false;
}

92 93 94 95 96 97 98 99 100 101
bool Monitor::applyToReferenceParticle(const Vector_t &R,
                                       const Vector_t &P,
                                       const double &t,
                                       Vector_t &,
                                       Vector_t &) {
    if (!OpalData::getInstance()->isInPrepState()) {
        const double dt = RefPartBunch_m->getdT();
        const double recpgamma = Physics::c * dt / Util::getGamma(P);
        const double middle = 0.5 * getElementLength();

102 103
        if (dt * R(2) < dt * middle && // multiply with dt to allow back tracking
            dt * (R(2) + P(2) * recpgamma) > dt * middle) {
104 105
            double frac = (middle - R(2)) / (P(2) * recpgamma);
            double time = t + frac * dt;
106
            Vector_t dR = frac * P * recpgamma;
Christof Metzger-Kraus's avatar
Christof Metzger-Kraus committed
107
            double ds = euclidean_norm(dR);
108 109 110 111 112 113 114 115 116 117 118
            lossDs_m->addReferenceParticle(csTrafoGlobal2Local_m.transformFrom(R + dR),
                                           csTrafoGlobal2Local_m.rotateFrom(P),
                                           time,
                                           RefPartBunch_m->get_sPos() + ds,
                                           RefPartBunch_m->getGlobalTrackStep());

            if (type_m == TEMPORAL) {
                const unsigned int localNum = RefPartBunch_m->getLocalNum();

                for (unsigned int i = 0; i < localNum; ++ i) {
                    const double recpgamma = Physics::c * dt / Util::getGamma(RefPartBunch_m->P[i]);
119 120 121 122 123 124
                    Vector_t shift = frac * recpgamma * RefPartBunch_m->P[i] - Vector_t(0, 0, middle);
                    lossDs_m->addParticle(RefPartBunch_m->R[i] + shift,
                                          RefPartBunch_m->P[i],
                                          RefPartBunch_m->ID[i],
                                          time,
                                          0);
125
                }
126
                OpalData::OPENMODE openMode;
127
                if (numPassages_m > 0) {
128 129 130
                    openMode = OpalData::OPENMODE::APPEND;
                } else {
                    openMode = OpalData::getInstance()->getOpenMode();
131
                }
132
                lossDs_m->save(1, openMode);
133 134 135 136 137
            }

            ++ numPassages_m;
        }
    }
gsell's avatar
gsell committed
138 139 140
    return false;
}

141
void Monitor::initialise(PartBunchBase<double, 3> *bunch, double &startField, double &endField) {
gsell's avatar
gsell committed
142
    RefPartBunch_m = bunch;
143 144 145
    endField = startField + halfLength_s;
    startField -= halfLength_s;

adelmann's avatar
adelmann committed
146
    if (filename_m == std::string(""))
147
        filename_m = getName();
adelmann's avatar
adelmann committed
148
    else
149 150 151 152 153 154 155 156
        filename_m = filename_m.substr(0, filename_m.rfind("."));

    const size_t totalNum = bunch->getTotalNum();
    double currentPosition = endField;
    if (totalNum > 0) {
        currentPosition = bunch->get_sPos();
    }

157 158
    if (OpalData::getInstance()->getOpenMode() == OpalData::OPENMODE::WRITE ||
        currentPosition < startField) {
159 160 161 162 163 164 165
        namespace fs = boost::filesystem;

        fs::path lossFileName = fs::path(filename_m + ".h5");
        if (fs::exists(lossFileName)) {
            Ippl::Comm->barrier();
            if (Ippl::myNode() == 0)
                fs::remove(lossFileName);
adelmann's avatar
adelmann committed
166

167 168 169 170 171
            Ippl::Comm->barrier();
        }
    }

    lossDs_m = std::unique_ptr<LossDataSink>(new LossDataSink(filename_m, !Options::asciidump, getType()));
gsell's avatar
gsell committed
172 173 174 175 176 177
}

void Monitor::finalise() {

}

kraus's avatar
kraus committed
178
void Monitor::goOnline(const double &) {
gsell's avatar
gsell committed
179 180 181 182
    online_m = true;
}

void Monitor::goOffline() {
183 184 185 186 187
    auto stats = lossDs_m->computeStatistics(numPassages_m);
    for (auto &stat: stats) {
        statFileEntries_sm.insert(std::make_pair(stat.spos_m, stat));
    }

188 189 190
    if (type_m != TEMPORAL) {
        lossDs_m->save(numPassages_m);
    }
adelmann's avatar
adelmann committed
191 192 193 194 195 196
}

bool Monitor::bends() const {
    return false;
}

197
void Monitor::setOutputFN(std::string fn) {
adelmann's avatar
adelmann committed
198 199 200 201
    filename_m = fn;
}

void Monitor::getDimensions(double &zBegin, double &zEnd) const {
202 203
    zBegin = -halfLength_s;
    zEnd = halfLength_s;
adelmann's avatar
adelmann committed
204 205 206
}


207 208
ElementBase::ElementType Monitor::getType() const {
    return MONITOR;
adelmann's avatar
adelmann committed
209 210
}

211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229
void Monitor::writeStatistics() {
    if (statFileEntries_sm.size() == 0) return;

    std::string fileName = OpalData::getInstance()->getInputBasename() + std::string("_Monitors.stat");
    auto instance = OpalData::getInstance();
    bool hasPriorTrack = instance->hasPriorTrack();
    bool inRestartRun = instance->inRestartRun();

    auto it = statFileEntries_sm.begin();
    double spos = it->first;
    Util::rewindLinesSDDS(fileName, spos, false);

    MonitorStatisticsWriter writer(fileName, hasPriorTrack || inRestartRun);

    for (const auto &entry: statFileEntries_sm) {
        writer.addRow(entry.second);
    }

    statFileEntries_sm.clear();
230
}