/Users/marc/Developer/EGSnrc/HEN_HOUSE/user_codes/egs_brachy/egs_brachy/muen.h

Go to the documentation of this file.

/*
################################################################################
#
#  egs_brachy muen.h
#  Copyright (C) 2016 Rowan Thomson, Dave Rogers, Randle Taylor, and Marc
#  Chamberland
#
#  This file is part of egs_brachy
#
#  egs_brachy is free software: you can redistribute it and/or modify it
#  under the terms of the GNU Affero General Public License as published
#  by the Free Software Foundation, either version 3 of the License, or
#  (at your option) any later version.
#
#  egs_brachy is distributed in the hope that it will be useful, but
#  WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#  Affero General Public License for more details:
#  <http://www.gnu.org/licenses/>.
#
################################################################################
#
#  When egs_brachy is used for publications, please cite our paper:
#  M. J. P. Chamberland, R. E. P. Taylor, D. W. O. Rogers, and R. M. Thomson,
#  egs brachy: a versatile and fast Monte Carlo code for brachytherapy,
#  Phys. Med. Biol. 61, 8214-8231 (2016).
#
################################################################################
#
#  Author:        Randle Taylor, 2016
#
#  Contributors:  Marc Chamberland
#                 Dave Rogers
#                 Rowan Thomson
#
################################################################################
*/
 
#include <algorithm>
#include <fstream>
#include <string>
#include <istream>
#include <iostream>
#include <sstream>
#include <map>
#include <cstdlib>
#include <vector>
 
#include "egs_functions.h"
#include "egs_input.h"
#include "egs_interpolator.h"
 
 
#ifndef MUENDATA_
#define MUENDATA_
 
namespace muen {
 
using namespace std;
 
 
std::vector<std::string> &split(const std::string &s, char delim, std::vector<std::string> &elems) {
    std::stringstream ss(s);
    std::string item;
    while (std::getline(ss, item, delim)) {
        elems.push_back(item);
    }
    return elems;
}
 
 
std::vector<std::string> split(const std::string &s, char delim) {
    std::vector<std::string> elems;
    split(s, delim, elems);
    return elems;
}
 
static inline std::string &ltrim(std::string &s) {
    s.erase(s.begin(), std::find_if(s.begin(), s.end(), [](unsigned char ch) {return !std::isspace(ch);}));
    return s;
}
 
static inline std::string &rtrim(std::string &s) {
    s.erase(std::find_if(s.rbegin(), s.rend(), [](unsigned char ch) {return !std::isspace(ch);}).base(), s.end());
    return s;
}
 
static inline std::string &trim(std::string &s) {
    return ltrim(rtrim(s));
}
 
 
typedef pair<double, double> MuenAtET;
 
typedef map<string, vector<MuenAtET> > MuenMapT;
 
class MuenDataParser {
 
    MuenMapT med_data; 
    static const int NSKIP = 3; 
    MuenMapT splitFileByMed(ifstream &in) {
 
        MuenMapT med_lines;
 
        string cur_med_name, dummy;
 
        while (in) {
 
            string line;
            getline(in, line);
            if (trim(line).size() == 0) {
                continue;
            }
 
            bool new_med = line.find(MUEN_START) != string::npos;
 
            if (new_med) {
                /* new medium found. Update the current medium name
                 * and move file pointer to start of data */
                cur_med_name = trim(split(line, '=')[1]);
                for (int i=0; i < NSKIP; i++) {
                    getline(in, dummy);
                    if (trim(dummy).size() == 0) {
                        i--;
                    }
                }
            } else {
 
                istringstream iss(line);
 
                double energy, muen;
                iss >> energy;
                iss >> dummy; // skip the MeV unit
                iss >> muen;
 
                if (iss.fail() && in) {
                    egsWarning("MuenDataParser:: Invalid muen data file format. Last line was:\n%s\n ", line.c_str());
                    med_lines.clear();
                    return med_lines;
                } else {
                    med_lines[cur_med_name].push_back(MuenAtET(energy, muen));
                }
            }
        }
 
        return med_lines;
 
    }
 
 
public:
 
    static const string MUEN_START; 
    MuenDataParser() {};
 
    int setMuenFile(string filename) {
 
        ifstream muen_data(filename.c_str());
 
        if (!muen_data) {
            return 1;
        }
 
        med_data = splitFileByMed(muen_data);
        muen_data.close();
        return 0;
 
    };
 
    EGS_Interpolator *getMuenInterpolator(string med_name) {
 
        size_t ndat = med_data[med_name].size();
        if (ndat < 2) {
            return 0;
        }
        EGS_Float emin = med_data[med_name][0].first;
        EGS_Float emax= med_data[med_name][ndat-1].first;
 
        EGS_Float *fmuen = new EGS_Float[ndat];
        for (size_t i=0; i < ndat; i++) {
            fmuen[i] = med_data[med_name][i].second;
        }
 
        EGS_Interpolator *muen = new EGS_Interpolator(ndat, log(emin), log(emax), fmuen);
        delete [] fmuen;
 
        return  muen;
 
    }
};
 
 
const string MuenDataParser::MUEN_START = "Muen values for medium MEDIUM =";
 
}
 
#endif