ProbeHistReader.cpp

Go to the documentation of this file.

//
// Class ProbeHistReader
//   Implements a parser and value extractor for hist files (*.hist).
//   It is for example used together with the septum objective.
//   A histogram file is generated by the OPAL probe element.
//
// Copyright (c) 2019, Matthias Frey, Paul Scherrer Institut, Villigen PSI, Switzerland
// All rights reserved
//
// Implemented as part of the PhD thesis
// "Precise Simulations of Multibunches in High Intensity Cyclotrons"
//
// 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 <iterator>
#include <regex>

#include "ProbeHistReader.h"
#include "Util/OptPilotException.h"

ProbeHistReader::ProbeHistReader(std::string filename)
    : filename_m(filename)
    , rmin_m(0.0)
    , binwidth_m(0.0)
    , bincount_m(0)
{ }


void ProbeHistReader::parseFile() {
    
    std::ifstream histfile;
    
    histfile.open(filename_m.c_str(), std::ios::in);
    
    if (!histfile) {
        throw OptPilotException("ProbeHistReader::parseFile()",
                                "Error opening file " + filename_m);
    }

    parseHeader(histfile);

    std::istream_iterator<size_t> it(histfile);
    while ( it != std::istream_iterator<size_t>() ) {
        bincount_m.push_back(*it);
        ++it;
    }

    histfile.close();

    if (histfile.is_open()) {
        throw OptPilotException("ProbeHistReader::parseFile()",
                                "Error closing file " + filename_m);
    }
}


size_t ProbeHistReader::minimum(double lower, double upper) {
    size_t lidx = (lower - rmin_m) / binwidth_m;
    size_t uidx = (upper - rmin_m) / binwidth_m;

    if (lidx >= uidx) {
        throw OptPilotException("ProbeHistReader::minimum()",
                                "Lower index >= upper index: " + std::to_string(lidx) +
                                " >= " + std::to_string(uidx));
    }

    if (uidx >= bincount_m.size()) {
        throw OptPilotException("ProbeHistReader::minimum()",
                                "Index >= number of bins: " + std::to_string(uidx) +
                                " >= " + std::to_string(bincount_m.size()));
    }

    container_t::iterator beg = std::begin(bincount_m);
    std::advance(beg, lidx);

    container_t::iterator end = std::begin(bincount_m);
    std::advance(end, uidx);

    container_t::iterator result = std::min_element(beg, end);

    if (result == bincount_m.end()) {
        throw OptPilotException("ProbeHistReader::minimum()",
                                "No minimum between " + std::to_string(lower) +
                                " and " + std::to_string(upper) + " found.");
    }

    return double(*result);
}


void ProbeHistReader::parseHeader(std::ifstream& ifs) {
    std::string header;
    std::getline(ifs, header);

    if (header.find("# Histogram bin counts") == std::string::npos) {
        throw OptPilotException("ProbeHistReader::parseHeader()",
                                "Error reading file " + filename_m);
    }

    const std::regex re("\\.*\\) (.*) mm (.*) mm (.*) (.*) mm");
    std::smatch match;
    std::regex_search(header, match, re);

    if ( match.size() != 5 ) {
        throw OptPilotException("ProbeHistReader::parseHeader()",
                                "Error parsing header of file " + filename_m);
    }

    rmin_m       = getValue<double>(match[1].str());
    double rmax  = getValue<double>(match[2].str());
    size_t nbins = getValue<size_t>(match[3].str());
    binwidth_m   = getValue<double>(match[4].str());

    if ( rmin_m >= rmax ) {
        throw OptPilotException("ProbeHistReader::parseHeader()",
                                "Not a valid histogram file " + filename_m);
    }

    bincount_m.reserve(nbins);
}