LossDataSink.cpp

Go to the documentation of this file.

//
// Class LossDataSink
//   This class writes file attributes to describe phase space of loss files
//
// Copyright (c) 200x - 2020, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// This file is part of OPAL.
//
// OPAL is free software: you can redistribute it and/or modify
// it under the terms of the GNU General Public License as published by
// the Free Software Foundation, either version 3 of the License, or
// (at your option) any later version.
//
// You should have received a copy of the GNU General Public License
// along with OPAL. If not, see <https://www.gnu.org/licenses/>.
//
#include "Structure/LossDataSink.h"
 
#include "AbstractObjects/OpalData.h"
#include "Algorithms/DistributionMoments.h"
#include "Message/GlobalComm.h"
#include "OPALconfig.h"
#include "Utilities/GeneralClassicException.h"
#include "Utilities/Options.h"
#include "Utilities/Util.h"
#include "Utility/IpplInfo.h"
 
#include <boost/filesystem.hpp>
 
#include <cmath>
 
extern Inform* gmsg;
 
#define WRITE_FILEATTRIB_STRING( attribute, value ) {             \
    h5_int64_t h5err = H5WriteFileAttribString (H5file_m, attribute, value); \
    if (h5err <= H5_ERR) {                                              \
        std::stringstream ss;                                               \
        ss << "failed to write string attribute " << attribute << " to file " << fn_m; \
        throw GeneralClassicException(std::string(__func__), ss.str()); \
    }\
}
#define WRITE_STEPATTRIB_FLOAT64( attribute, value, size ) { \
        h5_int64_t h5err = H5WriteStepAttribFloat64 (H5file_m, attribute, value, size); \
    if (h5err <= H5_ERR) { \
        std::stringstream ss; \
        ss << "failed to write float64 attribute " << attribute << " to file " << fn_m; \
        throw GeneralClassicException(std::string(__func__), ss.str()); \
    }\
}
#define WRITE_STEPATTRIB_INT64( attribute, value, size ) { \
        h5_int64_t h5err = H5WriteStepAttribInt64 (H5file_m, attribute, value, size); \
    if (h5err <= H5_ERR) { \
        std::stringstream ss; \
        ss << "failed to write int64 attribute " << attribute << " to file " << fn_m; \
        throw GeneralClassicException(std::string(__func__), ss.str()); \
    }\
}
#define WRITE_DATA_FLOAT64( name, value ) {                             \
        h5_int64_t h5err = H5PartWriteDataFloat64 (H5file_m, name, value); \
    if (h5err <= H5_ERR) { \
        std::stringstream ss; \
        ss << "failed to write float64 data " << name << " to file " << fn_m; \
        throw GeneralClassicException(std::string(__func__), ss.str()); \
    }\
}
#define WRITE_DATA_INT64( name, value ) { \
        h5_int64_t h5err = H5PartWriteDataInt64 (H5file_m, name, value); \
    if (h5err <= H5_ERR) { \
        std::stringstream ss; \
        ss << "failed to write int64 data " << name << " to file " << fn_m; \
        throw GeneralClassicException(std::string(__func__), ss.str()); \
    }\
}
#define SET_STEP() {               \
    h5_int64_t h5err = H5SetStep (H5file_m, H5call_m); \
    if (h5err <= H5_ERR) {                                              \
        std::stringstream ss;                                               \
        ss << "failed to set step " << H5call_m << " in file " << fn_m; \
        throw GeneralClassicException(std::string(__func__), ss.str()); \
    }\
}
#define GET_NUM_STEPS() {                             \
    H5call_m = H5GetNumSteps( H5file_m );                        \
    if (H5call_m <= H5_ERR) {                                             \
        std::stringstream ss;                                               \
        ss << "failed to get number of steps of file " << fn_m; \
        throw GeneralClassicException(std::string(__func__), ss.str()); \
    }\
}
#define SET_NUM_PARTICLES( num ) {             \
    h5_int64_t h5err = H5PartSetNumParticles (H5file_m, num); \
    if (h5err <= H5_ERR) {                                              \
        std::stringstream ss;                                               \
        ss << "failed to set number of particles to " << num << " in step " << H5call_m << " in file " << fn_m; \
        throw GeneralClassicException(std::string(__func__), ss.str()); \
    }\
}
 
#define OPEN_FILE( fname, mode, props ) {        \
    H5file_m = H5OpenFile (fname, mode, props); \
    if (H5file_m == (h5_file_t)H5_ERR) { \
        std::stringstream ss;                                               \
        ss << "failed to open file " << fn_m; \
        throw GeneralClassicException(std::string(__func__), ss.str()); \
    }\
}
#define CLOSE_FILE() {               \
    h5_int64_t h5err = H5CloseFile (H5file_m); \
    if (h5err <= H5_ERR) {                                              \
        std::stringstream ss;                                               \
        ss << "failed to close file " << fn_m; \
        throw GeneralClassicException(std::string(__func__), ss.str()); \
    }\
}
 
namespace {
    void f64transform(const std::vector<OpalParticle> &particles,
                   unsigned int startIdx,
                   unsigned int numParticles,
                   h5_float64_t *buffer,
                   std::function<h5_float64_t(const OpalParticle&)> select) {
        std::transform(particles.begin() + startIdx,
                       particles.begin() + startIdx + numParticles,
                       buffer, select);
    }
    void i64transform(const std::vector<OpalParticle> &particles,
                   unsigned int startIdx,
                   unsigned int numParticles,
                   h5_int64_t *buffer,
                   std::function<h5_int64_t(const OpalParticle&)> select) {
        std::transform(particles.begin() + startIdx,
                       particles.begin() + startIdx + numParticles,
                       buffer, select);
    }
 
    void cminmax(double &min, double &max, double val) {
        if (-val > min) {
            min = -val;
        } else if (val > max) {
            max = val;
        }
    }
}
 
SetStatistics::SetStatistics():
    outputName_m(""),
    spos_m(0.0),
    refTime_m(0.0),
    tmean_m(0.0),
    trms_m(0.0),
    nTotal_m(0),
    RefPartR_m(0.0),
    RefPartP_m(0.0),
    rmin_m(0.0),
    rmax_m(0.0),
    rmean_m(0.0),
    pmean_m(0.0),
    rrms_m(0.0),
    prms_m(0.0),
    rprms_m(0.0),
    normEmit_m(0.0),
    rsqsum_m(0.0),
    psqsum_m(0.0),
    rpsum_m(0.0),
    eps2_m(0.0),
    eps_norm_m(0.0),
    fac_m(0.0)
    { }
 
LossDataSink::LossDataSink(std::string outfn, bool hdf5Save, CollectionType collectionType):
    h5hut_mode_m(hdf5Save),
    H5file_m(0),
    outputName_m(outfn),
    H5call_m(0),
    collectionType_m(collectionType)
{
    particles_m.clear();
    turnNumber_m.clear();
    bunchNumber_m.clear();
 
    if (h5hut_mode_m && !Options::enableHDF5) {
        throw GeneralClassicException(
                "LossDataSink::LossDataSink",
                "You must select an OPTION to save Loss data files\n"
                "Please, choose 'ENABLEHDF5=TRUE' or 'ASCIIDUMP=TRUE'");
    }
 
    OpalData::getInstance()->checkAndAddOutputFileName(outputName_m);
 
    if (OpalData::getInstance()->hasBunchAllocated()) {
        OpalData::getInstance()->setOpenMode(OpalData::OPENMODE::APPEND);
    }
}
 
LossDataSink::LossDataSink(const LossDataSink &rhs):
    h5hut_mode_m(rhs.h5hut_mode_m),
    H5file_m(rhs.H5file_m),
    outputName_m(rhs.outputName_m),
    H5call_m(rhs.H5call_m),
    RefPartR_m(rhs.RefPartR_m),
    RefPartP_m(rhs.RefPartP_m),
    globalTrackStep_m(rhs.globalTrackStep_m),
    refTime_m(rhs.refTime_m),
    spos_m(rhs.spos_m),
    collectionType_m(rhs.collectionType_m)
{
    particles_m.clear();
    turnNumber_m.clear();
    bunchNumber_m.clear();
}
 
LossDataSink::~LossDataSink() noexcept(false) {
    if (H5file_m) {
        CLOSE_FILE ();
        H5file_m = 0;
    }
    Ippl::Comm->barrier();
}
 
void LossDataSink::openH5(h5_int32_t mode) {
    h5_prop_t props = H5CreateFileProp ();
    MPI_Comm comm = Ippl::getComm();
    H5SetPropFileMPIOCollective (props, &comm);
    OPEN_FILE (fn_m.c_str(), mode, props);
    H5CloseProp (props);
}
 
 
void LossDataSink::writeHeaderH5() {
 
    // Write file attributes to describe phase space to H5 file.
    std::stringstream OPAL_version;
    OPAL_version << OPAL_PROJECT_NAME << " " << OPAL_PROJECT_VERSION << " # git rev. " << Util::getGitRevision();
    WRITE_FILEATTRIB_STRING("OPAL_version", OPAL_version.str().c_str());
 
    WRITE_FILEATTRIB_STRING("SPOSUnit", "m");
    WRITE_FILEATTRIB_STRING("TIMEUnit", "s");
    WRITE_FILEATTRIB_STRING("RefPartRUnit", "m");
    WRITE_FILEATTRIB_STRING("RefPartPUnit", "#beta#gamma");
    WRITE_FILEATTRIB_STRING("GlobalTrackStepUnit", "1");
 
    WRITE_FILEATTRIB_STRING("centroidUnit", "m");
    WRITE_FILEATTRIB_STRING("RMSXUnit", "m");
    WRITE_FILEATTRIB_STRING("MEANPUnit", "#beta#gamma");
    WRITE_FILEATTRIB_STRING("RMSPUnit", "#beta#gamma");
    WRITE_FILEATTRIB_STRING("#varepsilonUnit", "m rad");
    WRITE_FILEATTRIB_STRING("#varepsilon-geomUnit", "m rad");
    WRITE_FILEATTRIB_STRING("ENERGYUnit", "MeV");
    WRITE_FILEATTRIB_STRING("dEUnit", "MeV");
    WRITE_FILEATTRIB_STRING("TotalChargeUnit", "C");
    WRITE_FILEATTRIB_STRING("TotalMassUnit", "MeV");
 
    WRITE_FILEATTRIB_STRING("idUnit", "1");
    WRITE_FILEATTRIB_STRING("xUnit", "m");
    WRITE_FILEATTRIB_STRING("yUnit", "m");
    WRITE_FILEATTRIB_STRING("zUnit", "m");
    WRITE_FILEATTRIB_STRING("pxUnit", "#beta#gamma");
    WRITE_FILEATTRIB_STRING("pyUnit", "#beta#gamma");
    WRITE_FILEATTRIB_STRING("pzUnit", "#beta#gamma");
    WRITE_FILEATTRIB_STRING("qUnit", "C");
    WRITE_FILEATTRIB_STRING("mUnit", "MeV");
 
    WRITE_FILEATTRIB_STRING("turnUnit", "1");
    WRITE_FILEATTRIB_STRING("bunchNumUnit", "1");
 
    WRITE_FILEATTRIB_STRING("timeUnit", "s");
    WRITE_FILEATTRIB_STRING("meanTimeUnit", "s");
    WRITE_FILEATTRIB_STRING("rmsTimeUnit", "s");
 
    if (collectionType_m == CollectionType::TEMPORAL) {
        WRITE_FILEATTRIB_STRING("type", "temporal");
    } else {
        WRITE_FILEATTRIB_STRING("type", "spatial");
    }
}
 
 
void LossDataSink::writeHeaderASCII() {
    bool hasTurn = hasTurnInformations();
    if (Ippl::myNode() == 0) {
        os_m << "# x (m),  y (m),  z (m),  px ( ),  py ( ),  pz ( ), id";
        if (hasTurn) {
            os_m << ",  turn ( ), bunchNumber ( )";
        }
        os_m << ", time (s)" << std::endl;
    }
}
 
void LossDataSink::addReferenceParticle(const Vector_t &x,
                                        const  Vector_t &p,
                                        double time,
                                        double spos,
                                        long long globalTrackStep
                                        ) {
    RefPartR_m.push_back(x);
    RefPartP_m.push_back(p);
    globalTrackStep_m.push_back((h5_int64_t)globalTrackStep);
    spos_m.push_back(spos);
    refTime_m.push_back(time);
}
 
void LossDataSink::addParticle(const OpalParticle &particle, const boost::optional<std::pair<int, short>> &turnBunchNumPair) {
    if (turnBunchNumPair) {
        if (!particles_m.empty() && turnNumber_m.empty()) {
            throw GeneralClassicException("LossDataSink::addParticle",
                                          "Either no particle or all have turn number and bunch number");
        }
        turnNumber_m.push_back(turnBunchNumPair.get().first);
        bunchNumber_m.push_back(turnBunchNumPair.get().second);
    }
 
    particles_m.push_back(particle);
}
 
void LossDataSink::save(unsigned int numSets, OpalData::OPENMODE openMode) {
 
    if (outputName_m.empty()) return;
    if (hasNoParticlesToDump()) return;
 
    if (openMode == OpalData::OPENMODE::UNDEFINED) {
        openMode = OpalData::getInstance()->getOpenMode();
    }
 
    namespace fs = boost::filesystem;
    if (h5hut_mode_m) {
 
        fn_m = outputName_m + std::string(".h5");
        *gmsg << level2 << "Save " << fn_m << endl;
        if (openMode == OpalData::OPENMODE::WRITE || !fs::exists(fn_m)) {
            openH5();
            writeHeaderH5();
        } else {
            openH5(H5_O_APPENDONLY);
            GET_NUM_STEPS ();
        }
 
        for (unsigned int i = 0; i < numSets; ++ i) {
            saveH5(i);
        }
        CLOSE_FILE ();
        H5file_m = 0;
 
    } else {
 
        fn_m = outputName_m + std::string(".loss");
        *gmsg << level2 << "Save " << fn_m << endl;
        if (openMode == OpalData::OPENMODE::WRITE || !fs::exists(fn_m)) {
            openASCII();
            writeHeaderASCII();
        } else {
            appendASCII();
        }
        saveASCII();
        closeASCII();
    }
 
    Ippl::Comm->barrier();
 
    particles_m.clear();
    turnNumber_m.clear();
    bunchNumber_m.clear();
    spos_m.clear();
    refTime_m.clear();
    RefPartR_m.clear();
    RefPartP_m.clear();
    globalTrackStep_m.clear();
}
 
// Note: This was changed to calculate the global number of dumped particles
// because there are two cases to be considered:
// 1. ALL nodes have 0 lost particles -> nothing to be done.
// 2. Some nodes have 0 lost particles, some not -> H5 can handle that but all
// nodes HAVE to participate, otherwise H5 waits endlessly for a response from
// the nodes that didn't enter the saveH5 function. -DW
bool LossDataSink::hasNoParticlesToDump() const {
 
    size_t nLoc = particles_m.size();
 
    reduce(nLoc, nLoc, OpAddAssign());
 
    return nLoc == 0;
}
 
bool LossDataSink::hasTurnInformations() const {
 
    bool hasTurnInformation = !turnNumber_m.empty();
 
    allreduce(hasTurnInformation, 1, std::logical_or<bool>());
 
    return hasTurnInformation > 0;
}
 
void LossDataSink::saveH5(unsigned int setIdx) {
    size_t startIdx = 0, endIdx = 0;
    size_t nLoc = particles_m.size();
    if (setIdx + 1 < startSet_m.size()) {
        startIdx = startSet_m[setIdx];
        endIdx = startSet_m[setIdx + 1];
        nLoc = endIdx - startIdx;
    }
 
    std::unique_ptr<size_t[]>  locN(new size_t[Ippl::getNodes()]);
    std::unique_ptr<size_t[]> globN(new size_t[Ippl::getNodes()]);
 
    for (int i = 0; i < Ippl::getNodes(); i++) {
        globN[i] = locN[i] = 0;
    }
 
    locN[Ippl::myNode()] = nLoc;
    reduce(locN.get(), locN.get() + Ippl::getNodes(), globN.get(), OpAddAssign());
 
    DistributionMoments engine;
    engine.compute(particles_m.begin() + startIdx, particles_m.begin() + endIdx);
 
    SET_STEP ();
    SET_NUM_PARTICLES (nLoc);
 
    if (setIdx < spos_m.size()) {
        WRITE_STEPATTRIB_FLOAT64 ("SPOS", &(spos_m[setIdx]), 1);
        WRITE_STEPATTRIB_FLOAT64 ("TIME", &(refTime_m[setIdx]), 1);
        WRITE_STEPATTRIB_FLOAT64 ("RefPartR", (h5_float64_t *)&(RefPartR_m[setIdx]), 3);
        WRITE_STEPATTRIB_FLOAT64 ("RefPartP", (h5_float64_t *)&(RefPartP_m[setIdx]), 3);
        WRITE_STEPATTRIB_INT64("GlobalTrackStep", &(globalTrackStep_m[setIdx]), 1);
    }
 
    Vector_t tmpVector;
    double tmpDouble;
    WRITE_STEPATTRIB_FLOAT64("centroid", (tmpVector = engine.getMeanPosition(), &tmpVector[0]), 3);
    WRITE_STEPATTRIB_FLOAT64("RMSX", (tmpVector = engine.getStandardDeviationPosition(), &tmpVector[0]), 3);
    WRITE_STEPATTRIB_FLOAT64("MEANP", (tmpVector = engine.getMeanMomentum(), &tmpVector[0]), 3);
    WRITE_STEPATTRIB_FLOAT64("RMSP", (tmpVector = engine.getStandardDeviationMomentum(), &tmpVector[0]), 3);
    WRITE_STEPATTRIB_FLOAT64("#varepsilon", (tmpVector = engine.getNormalizedEmittance(), &tmpVector[0]), 3);
    WRITE_STEPATTRIB_FLOAT64("#varepsilon-geom", (tmpVector = engine.getGeometricEmittance(), &tmpVector[0]), 3);
    WRITE_STEPATTRIB_FLOAT64("ENERGY", (tmpDouble = engine.getMeanKineticEnergy(), &tmpDouble), 1);
    WRITE_STEPATTRIB_FLOAT64("dE", (tmpDouble = engine.getStdKineticEnergy(), &tmpDouble), 1);
    WRITE_STEPATTRIB_FLOAT64("TotalCharge", (tmpDouble = engine.getTotalCharge(), &tmpDouble), 1);
    WRITE_STEPATTRIB_FLOAT64("TotalMass", (tmpDouble = engine.getTotalMass(), &tmpDouble), 1);
    WRITE_STEPATTRIB_FLOAT64("meanTime", (tmpDouble = engine.getMeanTime(), &tmpDouble), 1);
    WRITE_STEPATTRIB_FLOAT64("rmsTime", (tmpDouble = engine.getStdTime(), &tmpDouble), 1);
 
    // Write all data
    std::vector<char> buffer(nLoc * sizeof(h5_float64_t));
    h5_float64_t *f64buffer = nLoc > 0? reinterpret_cast<h5_float64_t*>(&buffer[0]) : nullptr;
    h5_int64_t *i64buffer = nLoc > 0? reinterpret_cast<h5_int64_t*>(&buffer[0]) : nullptr;
 
    ::i64transform(particles_m, startIdx, nLoc, i64buffer,
                   [](const OpalParticle &particle) { return particle.getId(); });
    WRITE_DATA_INT64 ("id", i64buffer);
    ::f64transform(particles_m, startIdx, nLoc, f64buffer,
                   [](const OpalParticle &particle) { return particle.getX(); });
    WRITE_DATA_FLOAT64 ("x", f64buffer);
    ::f64transform(particles_m, startIdx, nLoc, f64buffer,
                   [](const OpalParticle &particle) { return particle.getY(); });
    WRITE_DATA_FLOAT64 ("y", f64buffer);
    ::f64transform(particles_m, startIdx, nLoc, f64buffer,
                   [](const OpalParticle &particle) { return particle.getZ(); });
    WRITE_DATA_FLOAT64 ("z", f64buffer);
    ::f64transform(particles_m, startIdx, nLoc, f64buffer,
                   [](const OpalParticle &particle) { return particle.getPx(); });
    WRITE_DATA_FLOAT64 ("px", f64buffer);
    ::f64transform(particles_m, startIdx, nLoc, f64buffer,
                   [](const OpalParticle &particle) { return particle.getPy(); });
    WRITE_DATA_FLOAT64 ("py", f64buffer);
    ::f64transform(particles_m, startIdx, nLoc, f64buffer,
                   [](const OpalParticle &particle) { return particle.getPz(); });
    WRITE_DATA_FLOAT64 ("pz", f64buffer);
    ::f64transform(particles_m, startIdx, nLoc, f64buffer,
                   [](const OpalParticle &particle) { return particle.getCharge(); });
    WRITE_DATA_FLOAT64 ("q", f64buffer);
    ::f64transform(particles_m, startIdx, nLoc, f64buffer,
                   [](const OpalParticle &particle) { return particle.getMass(); });
    WRITE_DATA_FLOAT64 ("m", f64buffer);
 
    if (hasTurnInformations()) {
        std::copy(turnNumber_m.begin() + startIdx,
                  turnNumber_m.begin() + startIdx + nLoc,
                  i64buffer);
        WRITE_DATA_INT64 ("turn", i64buffer);
 
        std::copy(bunchNumber_m.begin() + startIdx,
                  bunchNumber_m.begin() + startIdx + nLoc,
                  i64buffer);
        WRITE_DATA_INT64 ("bunchNumber", i64buffer);
    }
 
    ::f64transform(particles_m, startIdx, nLoc, f64buffer,
                   [](const OpalParticle &particle) { return particle.getTime(); });
    WRITE_DATA_FLOAT64 ("time", f64buffer);
 
    ++ H5call_m;
}
 
void LossDataSink::saveASCII() {
 
    /*
      ASCII output
    */
    int tag = Ippl::Comm->next_tag(IPPL_APP_TAG3, IPPL_APP_CYCLE);
    bool hasTurn = hasTurnInformations();
    if (Ippl::Comm->myNode() == 0) {
        const unsigned partCount = particles_m.size();
 
        for (unsigned i = 0; i < partCount; i++) {
            const OpalParticle& particle = particles_m[i];
            os_m << particle.getX() << "   ";
            os_m << particle.getY() << "   ";
            os_m << particle.getZ() << "   ";
            os_m << particle.getPx() << "   ";
            os_m << particle.getPy() << "   ";
            os_m << particle.getPz() << "   ";
            os_m << particle.getId() << "   ";
            if (hasTurn) {
                os_m << turnNumber_m[i] << "   ";
                os_m << bunchNumber_m[i] << "   ";
            }
            os_m << particle.getTime()
                 << std::endl;
        }
 
        int notReceived =  Ippl::getNodes() - 1;
        while (notReceived > 0) {
            unsigned dataBlocks = 0;
            int node = COMM_ANY_NODE;
            Message *rmsg =  Ippl::Comm->receive_block(node, tag);
            if (rmsg == 0) {
                ERRORMSG("LossDataSink: Could not receive from client nodes output." << endl);
            }
            notReceived--;
            rmsg->get(&dataBlocks);
            rmsg->get(&hasTurn);
            for (unsigned i = 0; i < dataBlocks; i++) {
                long id;
                size_t bunchNum, turn;
                double rx, ry, rz, px, py, pz, time;
 
                os_m << (rmsg->get(&rx), rx) << "   ";
                os_m << (rmsg->get(&ry), ry) << "   ";
                os_m << (rmsg->get(&rz), rz) << "   ";
                os_m << (rmsg->get(&px), px) << "   ";
                os_m << (rmsg->get(&py), py) << "   ";
                os_m << (rmsg->get(&pz), pz) << "   ";
                os_m << (rmsg->get(&id), id) << "   ";
                if (hasTurn) {
                    os_m << (rmsg->get(&turn), turn) << "   ";
                    os_m << (rmsg->get(&bunchNum), bunchNum) << "   ";
                }
                os_m << (rmsg->get(&time), time)
                     << std::endl;
            }
            delete rmsg;
        }
    } else {
        Message *smsg = new Message();
        const unsigned msgsize = particles_m.size();
        smsg->put(msgsize);
        smsg->put(hasTurn);
        for (unsigned i = 0; i < msgsize; i++) {
            const OpalParticle& particle = particles_m[i];
            smsg->put(particle.getX());
            smsg->put(particle.getY());
            smsg->put(particle.getZ());
            smsg->put(particle.getPx());
            smsg->put(particle.getPy());
            smsg->put(particle.getPz());
            smsg->put(particle.getId());
            if (hasTurn) {
                smsg->put(turnNumber_m[i]);
                smsg->put(bunchNumber_m[i]);
            }
            smsg->put(particle.getTime());
        }
        bool res = Ippl::Comm->send(smsg, 0, tag);
        if (!res) {
            ERRORMSG("LossDataSink Ippl::Comm->send(smsg, 0, tag) failed " << endl);
        }
    }
}
 
void LossDataSink::splitSets(unsigned int numSets) {
    if (numSets <= 1 ||
        particles_m.size() == 0) return;
 
    const size_t nLoc = particles_m.size();
    size_t avgNumPerSet = nLoc / numSets;
    std::vector<size_t> numPartsInSet(numSets, avgNumPerSet);
    for (unsigned int j = 0; j < (nLoc - numSets * avgNumPerSet); ++ j) {
        ++ numPartsInSet[j];
    }
    startSet_m.resize(numSets + 1, 0);
 
    std::vector<double> data(2 * numSets, 0.0);
    double* meanT = &data[0];
    double* numParticles = &data[numSets];
    std::vector<double> timeRange(numSets, 0.0);
    double maxT = particles_m[0].getTime();
 
    for (unsigned int iteration = 0; iteration < 2; ++ iteration) {
        size_t partIdx = 0;
        for (unsigned int j = 0; j < numSets; ++ j) {
 
            const size_t &numThisSet = numPartsInSet[j];
            for (size_t k = 0; k < numThisSet; ++ k, ++ partIdx) {
                meanT[j] += particles_m[partIdx].getTime();
                maxT = std::max(maxT, particles_m[partIdx].getTime());
            }
            numParticles[j] = numThisSet;
        }
 
        allreduce(&(data[0]), 2 * numSets, std::plus<double>());
 
        for (unsigned int j = 0; j < numSets; ++ j) {
            meanT[j] /= numParticles[j];
        }
 
        for (unsigned int j = 0; j < numSets - 1; ++ j) {
            timeRange[j] = 0.5 * (meanT[j] + meanT[j + 1]);
        }
        timeRange[numSets - 1] = maxT;
 
        std::fill(numPartsInSet.begin(),
                  numPartsInSet.end(),
                  0);
 
        size_t setNum = 0;
        size_t idxPrior = 0;
        for (size_t idx = 0; idx < nLoc; ++ idx) {
            if (particles_m[idx].getTime() > timeRange[setNum]) {
                numPartsInSet[setNum] = idx - idxPrior;
                idxPrior = idx;
                ++ setNum;
            }
        }
        numPartsInSet[numSets - 1] = nLoc - idxPrior;
    }
 
    for (unsigned int i = 0; i < numSets; ++ i) {
        startSet_m[i + 1] = startSet_m[i] + numPartsInSet[i];
    }
}
 
SetStatistics LossDataSink::computeSetStatistics(unsigned int setIdx) {
    SetStatistics stat;
    double part[6];
 
    const unsigned int totalSize = 45;
    double plainData[totalSize];
    double rminmax[6] = {};
 
    Util::KahanAccumulation data[totalSize];
    Util::KahanAccumulation *localCentroid = data + 1;
    Util::KahanAccumulation *localMoments = data + 7;
    Util::KahanAccumulation *localOthers = data + 43;
 
    size_t startIdx = 0;
    size_t nLoc = particles_m.size();
    if (setIdx + 1 < startSet_m.size()) {
        startIdx = startSet_m[setIdx];
        nLoc = startSet_m[setIdx + 1] - startSet_m[setIdx];
    }
 
    data[0].sum = nLoc;
 
    unsigned int idx = startIdx;
    for (unsigned long k = 0; k < nLoc; ++ k, ++ idx) {
        const OpalParticle& particle = particles_m[idx];
 
        part[0] = particle.getX();
        part[1] = particle.getPx();
        part[2] = particle.getY();
        part[3] = particle.getPy();
        part[4] = particle.getZ();
        part[5] = particle.getPz();
 
        for (int i = 0; i < 6; i++) {
            localCentroid[i] += part[i];
            for (int j = 0; j <= i; j++) {
                localMoments[i * 6 + j] += part[i] * part[j];
            }
        }
        localOthers[0] += particle.getTime();
        localOthers[1] += std::pow(particle.getTime(), 2);
 
        ::cminmax(rminmax[0], rminmax[1], particle.getX());
        ::cminmax(rminmax[2], rminmax[3], particle.getY());
        ::cminmax(rminmax[4], rminmax[5], particle.getZ());
    }
 
    for (int i = 0; i < 6; i++) {
        for (int j = 0; j < i; j++) {
            localMoments[j * 6 + i] = localMoments[i * 6 + j];
        }
    }
 
    for (unsigned int i = 0; i < totalSize; ++ i) {
        plainData[i] = data[i].sum;
    }
 
    new_reduce(plainData, totalSize, std::plus<double>());
    new_reduce(rminmax, 6, std::greater<double>());
 
    if (Ippl::myNode() != 0) return stat;
    if (plainData[0] == 0.0) return stat;
 
    double *centroid = plainData + 1;
    double *moments = plainData + 7;
    double *others = plainData + 43;
 
    stat.outputName_m = outputName_m;
    stat.spos_m = spos_m[setIdx];
    stat.refTime_m = refTime_m[setIdx];
    stat.RefPartR_m = RefPartR_m[setIdx];
    stat.RefPartP_m = RefPartP_m[setIdx];
    stat.nTotal_m = (unsigned long)std::round(plainData[0]);
 
    for (unsigned int i = 0 ; i < 3u; i++) {
        stat.rmean_m(i) = centroid[2 * i] / stat.nTotal_m;
        stat.pmean_m(i) = centroid[(2 * i) + 1] / stat.nTotal_m;
        stat.rsqsum_m(i) = (moments[2 * i * 6 + 2 * i] -
                            stat.nTotal_m * std::pow(stat.rmean_m(i), 2));
        stat.psqsum_m(i) = std::max(0.0,
                                    moments[(2 * i + 1) * 6 + (2 * i) + 1] -
                                    stat.nTotal_m * std::pow(stat.pmean_m(i), 2));
        stat.rpsum_m(i) = (moments[(2 * i) * 6 + (2 * i) + 1] -
                           stat.nTotal_m * stat.rmean_m(i) * stat.pmean_m(i));
    }
    stat.tmean_m = others[0] / stat.nTotal_m;
    stat.trms_m = std::sqrt(std::max(0.0, (others[1] / stat.nTotal_m - std::pow(stat.tmean_m, 2))));
 
    stat.eps2_m = ((stat.rsqsum_m * stat.psqsum_m - stat.rpsum_m * stat.rpsum_m) /
                   (1.0 * stat.nTotal_m * stat.nTotal_m));
 
    stat.rpsum_m /= stat.nTotal_m;
 
    for (unsigned int i = 0 ; i < 3u; i++) {
        stat.rrms_m(i) = std::sqrt(std::max(0.0, stat.rsqsum_m(i)) / stat.nTotal_m);
        stat.prms_m(i) = std::sqrt(std::max(0.0, stat.psqsum_m(i)) / stat.nTotal_m);
        stat.eps_norm_m(i)  =  std::sqrt(std::max(0.0, stat.eps2_m(i)));
        double tmp = stat.rrms_m(i) * stat.prms_m(i);
        stat.fac_m(i) = (tmp == 0) ? 0.0 : 1.0 / tmp;
        stat.rmin_m(i) = -rminmax[2 * i];
        stat.rmax_m(i) = rminmax[2 * i + 1];
    }
 
    stat.rprms_m = stat.rpsum_m * stat.fac_m;
 
    return stat;
}