Monitor.cpp 9.61 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/PartBunch.h"
gsell's avatar
gsell committed
23
#include "AbsBeamline/BeamlineVisitor.h"
adelmann's avatar
adelmann committed
24 25
#include "Fields/Fieldmap.hh"
#include "Structure/LossDataSink.h"
adelmann's avatar
adelmann committed
26 27
#include "Utilities/Options.h"

adelmann's avatar
adelmann committed
28
#include <memory>
29

adelmann's avatar
adelmann committed
30
using namespace Options;
31

gsell's avatar
gsell committed
32
#include <fstream>
33
#include <memory>
gsell's avatar
gsell committed
34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57

using namespace std;

// Class Monitor
// ------------------------------------------------------------------------

Monitor::Monitor():
    Component(),
    filename_m(""),
    plane_m(OFF),
    position_m(0.0),
    informed_m(false)
{}


Monitor::Monitor(const Monitor &right):
    Component(right),
    filename_m(right.filename_m),
    plane_m(right.plane_m),
    position_m(right.position_m),
    informed_m(right.informed_m)
{}


58
Monitor::Monitor(const std::string &name):
gsell's avatar
gsell committed
59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74
    Component(name),
    filename_m(""),
    plane_m(OFF),
    position_m(0.0),
    informed_m(false)
{}


Monitor::~Monitor()
{}


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

75
bool Monitor::apply(const size_t &i, const double &t, double E[], double B[]) {
gsell's avatar
gsell committed
76 77 78 79
    Vector_t Ev(0, 0, 0), Bv(0, 0, 0);
    return apply(i, t, Ev, Bv);
}

80
bool Monitor::apply(const size_t &i, const double &t, Vector_t &E, Vector_t &B) {
gsell's avatar
gsell committed
81 82 83 84 85
    const Vector_t &R = RefPartBunch_m->R[i];
    const Vector_t &P = RefPartBunch_m->P[i];
    const double recpgamma = Physics::c * RefPartBunch_m->getdT() / sqrt(1.0  + dot(P, P));
    if(online_m && R(2) < position_m && R(2) + P(2) * recpgamma > position_m) {
        double frac = (position_m - R(2)) / (P(2) * recpgamma);
adelmann's avatar
adelmann committed
86 87

        lossDs_m->addParticle(Vector_t(R(0) + frac * P(0) * recpgamma, R(1) + frac * P(1) * recpgamma, position_m),
kraus's avatar
kraus committed
88
                                   P, RefPartBunch_m->ID[i], t + frac * RefPartBunch_m->getdT(), 0);
gsell's avatar
gsell committed
89 90 91 92 93 94 95 96 97 98 99 100 101 102
    }

    return false;
}

bool Monitor::apply(const Vector_t &R, const Vector_t &centroid, const double &t, Vector_t &E, Vector_t &B) {
    return false;
}

void Monitor::initialise(PartBunch *bunch, double &startField, double &endField, const double &scaleFactor) {
    RefPartBunch_m = bunch;
    position_m = startField;
    startField -= 0.005;
    endField = position_m + 0.005;
adelmann's avatar
adelmann committed
103 104 105 106 107
    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));

gsell's avatar
gsell committed
108 109 110 111 112 113 114 115 116 117
}

void Monitor::finalise() {

}

void Monitor::goOnline() {
    if(RefPartBunch_m == NULL) {
        if(!informed_m) {
            Inform msg("Monitor ");
118
            std::string errormsg;
gsell's avatar
gsell committed
119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141
            errormsg = Fieldmap::typeset_msg("BUNCH SIZE NOT SET", "warning");
            msg << errormsg << "\n"
                << endl;
            if(Ippl::myNode() == 0) {
                ofstream omsg("errormsg.txt", ios_base::app);
                omsg << errormsg << endl;
                omsg.close();
            }
            informed_m = true;
        }
        return;
    }

    if(Monitor::h5pfiles_s.find(filename_m) == Monitor::h5pfiles_s.end()) {
        Monitor::h5pfiles_s.insert(pair<string, unsigned int>(filename_m, 1));
        step_m = 0;
    } else {
        step_m = (*Monitor::h5pfiles_s.find(filename_m)).second ++;
    }
    online_m = true;
}

void Monitor::goOffline() {
adelmann's avatar
adelmann committed
142 143 144 145 146 147 148
    lossDs_m->save();
}

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

149
void Monitor::setOutputFN(std::string fn) {
adelmann's avatar
adelmann committed
150 151 152 153 154 155 156 157 158
    filename_m = fn;
}

void Monitor::getDimensions(double &zBegin, double &zEnd) const {
    zBegin = position_m - 0.005;
    zEnd = position_m + 0.005;
}


159 160
const std::string &Monitor::getType() const {
    static const std::string type("Monitor");
adelmann's avatar
adelmann committed
161 162 163 164 165 166
    return type;
}

void Monitor::moveBy(const double &dz) {
    position_m += dz;
}
adelmann's avatar
adelmann committed
167

adelmann's avatar
adelmann committed
168 169 170 171 172 173
map<string, unsigned int> Monitor::h5pfiles_s = map<string, unsigned int>();




    /*
adelmann's avatar
adelmann committed
174
    if (!Options::enableHDF5) return;
kraus's avatar
kraus committed
175

adelmann's avatar
adelmann committed
176
	reduce(online_m, online_m, OpOr());
gsell's avatar
gsell committed
177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230

    if(online_m) {
        online_m = false;
        if(filename_m == "") return;

        unsigned long nLoc = PosX_m.size();
        unsigned long i = 0;
        h5_file_t *H5file;
        h5_int64_t rc;
        if(step_m == 0) {
#ifdef PARALLEL_IO
            H5file = H5OpenFile(filename_m.c_str(), H5_O_WRONLY, Ippl::getComm());
#else
            H5file = H5OpenFile(filename_m.c_str(), H5_O_WRONLY, 0);
#endif
            rc = H5WriteFileAttribString(H5file, "timeUnit", "s");
            if(rc != H5_SUCCESS)
                ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
            rc = H5WriteFileAttribString(H5file, "xUnit", "m");
            if(rc != H5_SUCCESS)
                ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
            rc = H5WriteFileAttribString(H5file, "yUnit", "m");
            if(rc != H5_SUCCESS)
                ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
            rc = H5WriteFileAttribString(H5file, "pxUnit", "#beta#gamma");
            if(rc != H5_SUCCESS)
                ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
            rc = H5WriteFileAttribString(H5file, "pyUnit", "#beta#gamma");
            if(rc != H5_SUCCESS)
                ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
            rc = H5WriteFileAttribString(H5file, "pzUnit", "#beta#gamma");
            if(rc != H5_SUCCESS)
                ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
            rc = H5WriteFileAttribString(H5file, "SPOSUnit", "m");
            if(rc != H5_SUCCESS)
                ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
        } else {
#ifdef PARALLEL_IO
            H5file = H5OpenFile(filename_m.c_str(), H5_O_APPEND, Ippl::getComm());
#else
            H5file = H5OpenFile(filename_m.c_str(), H5_O_APPEND, 0);
#endif
        }

        rc = H5SetStep(H5file, step_m);
        if(rc != H5_SUCCESS)
            ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
        rc = H5WriteStepAttribFloat64(H5file, "SPOS", &position_m, 1);
        if(rc != H5_SUCCESS)
            ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
        rc = H5PartSetNumParticles(H5file, PosX_m.size());
        if(rc != H5_SUCCESS)
            ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);

231 232 233
        std::unique_ptr<char> varray(new char[nLoc * sizeof(double)]);
        double *fvalues = reinterpret_cast<double*>(varray.get());
        h5_int64_t *ids = reinterpret_cast<h5_int64_t*>(varray.get());
kraus's avatar
kraus committed
234 235

	  FixMe: if I write with nLoc==0 -> rc == -2
adelmann's avatar
adelmann committed
236

237

gsell's avatar
gsell committed
238

239 240
	if (nLoc > 0) {
	  for(i = 0; i < nLoc; ++i) {
gsell's avatar
gsell committed
241 242
            fvalues[i] = PosX_m.front();
            PosX_m.pop_front();
243 244 245 246 247
	  }
	  rc = H5PartWriteDataFloat64(H5file, "x", fvalues);
	  if(rc != H5_SUCCESS)
	    ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << " nloc= " << nLoc << " fn= " << filename_m << endl);
	  for(i = 0; i < nLoc; ++i) {
gsell's avatar
gsell committed
248 249
            fvalues[i] = PosY_m.front();
            PosY_m.pop_front();
250 251 252
	  }
	  rc = H5PartWriteDataFloat64(H5file, "y", fvalues);
	  if(rc != H5_SUCCESS)
gsell's avatar
gsell committed
253
            ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
254
	  for(i = 0; i < nLoc; ++i) {
gsell's avatar
gsell committed
255 256
            fvalues[i] = MomentumX_m.front();
            MomentumX_m.pop_front();
257 258 259
	  }
	  rc = H5PartWriteDataFloat64(H5file, "px", fvalues);
	  if(rc != H5_SUCCESS)
gsell's avatar
gsell committed
260
            ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
261
	  for(i = 0; i < nLoc; ++i) {
gsell's avatar
gsell committed
262 263
            fvalues[i] = MomentumY_m.front();
            MomentumY_m.pop_front();
264 265 266
	  }
	  rc = H5PartWriteDataFloat64(H5file, "py", fvalues);
	  if(rc != H5_SUCCESS)
gsell's avatar
gsell committed
267
            ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
268
	  for(i = 0; i < nLoc; ++i) {
gsell's avatar
gsell committed
269 270
            fvalues[i] = MomentumZ_m.front();
            MomentumZ_m.pop_front();
271 272 273
	  }
	  rc = H5PartWriteDataFloat64(H5file, "pz", fvalues);
	  if(rc != H5_SUCCESS)
gsell's avatar
gsell committed
274
            ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
275
	  for(i = 0; i < nLoc; ++i) {
gsell's avatar
gsell committed
276 277
            fvalues[i] = time_m.front();
            time_m.pop_front();
278 279 280
	  }
	  rc = H5PartWriteDataFloat64(H5file, "time", fvalues);
	  if(rc != H5_SUCCESS)
gsell's avatar
gsell committed
281
            ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
282
	  for(i = 0; i < nLoc; ++i) {
gsell's avatar
gsell committed
283 284
            ids[i] = id_m.front();
            id_m.pop_front();
285 286 287
	  }
	  rc = H5PartWriteDataInt64(H5file, "id", ids);
	  if(rc != H5_SUCCESS)
gsell's avatar
gsell committed
288
            ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
289
	}
gsell's avatar
gsell committed
290 291
        rc = H5CloseFile(H5file);
        if(rc != H5_SUCCESS)
292
	  ERRORMSG("H5 rc= " << rc << " in " << __FILE__ << " @ line " << __LINE__ << endl);
gsell's avatar
gsell committed
293
    }
adelmann's avatar
adelmann committed
294
    */