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

Go to the documentation of this file.

/*
################################################################################
#
#  egs_brachy spec_scoring.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 <sstream>
#include "egs_advanced_application.h"
#include "egs_functions.h"
#include "egs_input.h"
#include "egs_interface2.h"
#include "egs_scoring.h"
#include "ginfo.h"
#include "pubsub.h"
 
 
class EnergyScoringStats : public Subscriber {
 
    EGS_Float total_energy_initialized; 
    EGS_Float energy_escaping_sources;  
    EGS_Float energy_escaping_geom;  
    EGS_Float getParticleEnergy(const EGS_Particle *p, bool subtractRM=true) {
        EGS_Float rest_mass = 0;
        if (subtractRM && p->q != 0){
            rest_mass = the_useful->prm;
        }
        return (p->E - rest_mass)*p->wt;
    }
 
    void scoreParticleInitialized(EGS_Particle *p) {
        /* When particle is first initialized by getNextParticle it does not
         * include rest mass so we tell getParticleEnergy not to subtract the
         * rest mass here. */
        total_energy_initialized += getParticleEnergy(p, false);
    }
 
    void scoreParticleEscapingSource(EGS_Particle *p) {
        energy_escaping_sources += getParticleEnergy(p);
    }
 
    void scoreParticleEscapingGeom(EGS_Particle *p) {
        energy_escaping_geom += getParticleEnergy(p);
    }
 
public:
 
    EnergyScoringStats(Publisher *publisher) :
        total_energy_initialized(0),
        energy_escaping_sources(0),
        energy_escaping_geom(0) {
 
        publisher->subscribe(this, PARTICLE_INITIALIZED);
        publisher->subscribe(this, PARTICLE_ESCAPING_SOURCE);
        publisher->subscribe(this, PARTICLE_ESCAPING_GEOM);
 
    };
 
 
    void scoreEnergyInitialized(EGS_Float E) {
        total_energy_initialized += E;
    }
    EGS_Float escapingSourcesRatio() {
        if (total_energy_initialized > 0) {
            return energy_escaping_sources/total_energy_initialized;
        }
        return 0;
    };
 
 
    EGS_Float escapingGeomRatio() {
        if (total_energy_initialized > 0) {
            return energy_escaping_geom/total_energy_initialized;
        }
 
        return 0;
    }
 
    EGS_Float totalEnergyInitialized() {
        return total_energy_initialized;
    }
 
    EGS_Float energyEscapingSources() {
        return energy_escaping_sources;
    }
 
    EGS_Float energyEscapingGeom() {
        return energy_escaping_geom;
    }
 
 
    void update(EB_Message message, void *particle) {
 
        if (message == PARTICLE_INITIALIZED) {
            EGS_Particle *p = static_cast<EGS_Particle *>(particle);
            scoreParticleInitialized(p);
        } else if (message == PARTICLE_ESCAPING_SOURCE) {
            EGS_Particle *p = static_cast<EGS_Particle *>(particle);
            scoreParticleEscapingSource(p);
        } else if (message == PARTICLE_ESCAPING_GEOM) {
 
            EGS_Particle *p = static_cast<EGS_Particle *>(particle);
            scoreParticleEscapingGeom(p);
        }
 
    }
 
 
    void outputResults() {
 
        string sep(80, '=');
 
        egsInformation("\n\nEnergy Initialized and Escaping\n%s\n", sep.c_str());
        egsInformation("Total Energy Initialized    : %.5G MeV", totalEnergyInitialized());
        egsInformation(
            "\nTotal Energy Escaping Sources : %.5G MeV (%.3f%%)",
            energyEscapingSources(),
            escapingSourcesRatio()*100
        );
 
        egsInformation(
            "\nTotal Energy Escaping Geometry: %.5G MeV (%.3f%%)\n",
            energyEscapingGeom(),
            escapingGeomRatio()*100
        );
 
    }
 
    int outputData(ostream *ofile) {
        (*ofile) << total_energy_initialized << " ";
        (*ofile) << energy_escaping_sources << " ";
        (*ofile) << energy_escaping_geom << "\n";
        return 0;
    };
 
    int readData(istream *ifile) {
        (*ifile) >> total_energy_initialized;
        (*ifile) >> energy_escaping_sources;
        (*ifile) >> energy_escaping_geom;
        return 0;
    }
 
    void resetCounter() {
        total_energy_initialized = 0;
        energy_escaping_sources = 0;
        energy_escaping_geom = 0;
    }
 
    int addState(istream &data) {
        EGS_Float total_energy_initialized_tmp;
        EGS_Float energy_escaping_sources_tmp;
        EGS_Float energy_escaping_geom_tmp;
 
        data >> total_energy_initialized_tmp;
        data >> energy_escaping_sources_tmp;
        data >> energy_escaping_geom_tmp;
 
        total_energy_initialized += total_energy_initialized_tmp;
        energy_escaping_sources += energy_escaping_sources_tmp;
        energy_escaping_geom += energy_escaping_geom_tmp;
 
        return 0;
    };
 
};
 
 
class BaseSpectrumScorer : public Subscriber {
 
protected:
 
    string format; 
    int egsnrc_mode; 
    int particle_type; 
    EGS_I64 cur_history; 
    EGS_I64 eff_history; 
    EGS_Float bin_width;
    EGS_Float total_scored;
 
    EGS_BaseSource *source; 
    bool valid;
 
    string fextension; 
    string outputCSV(string filename);
 
    string outputXMGR(string filename);
 
    string outputEGSnrc(string filename);
 
    int getBin(EGS_Float E) const;
 
    virtual string getTitle() const {
        return "EGS_Brachy: spectrum";
    }
 
    virtual string getSubTitle() const {
        return "";
    }
 
    virtual string getYAxisLabel() const {
        return "";
    }
 
    virtual string getXAxisLabel() const {
        return "E / MeV";
    }
 
    string getFileName(string root) const {
        return root + "." + getFileExtension();
    }
 
    virtual string getFileExtension() const {
        return "spec";
    };
 
    virtual void getResult(int bin, EGS_Float &r, EGS_Float &dr);
 
    string getParticleName() const {
 
        if (particle_type == 0) {
            return "photon";
        } else if (particle_type == -1) {
            return "electron";
        } else if (particle_type == 1) {
            return "positron";
        }
 
        return "";
    }
 
    virtual void outputTotal() {
        egsInformation("Total quantity scored: %.3g\n", total_scored);
    };
 
 
public:
 
    static const int DEFAULT_NBINS = 100;
 
 
    int nbins; 
    EGS_ScoringArray *bins; 
    EGS_Float e_min; 
    EGS_Float e_max; 
    BaseSpectrumScorer(EGS_Input *input, EGS_BaseSource *src, GeomInfo *ginfo, Publisher *publisher);
 
    static BaseSpectrumScorer *getSpectrumScorer(EGS_Input *inp, EGS_BaseSource *source, GeomInfo *ginfo,
            Publisher *publisher);
 
    virtual ~BaseSpectrumScorer();
 
    EGS_Float getBinWidth() const;
 
    double getParticleEnergy(const EGS_Particle *p) const;
 
    virtual void update(EB_Message message, void *data=0);
 
    virtual void score(EB_Message message, void *data=0) = 0;
 
    void setEffectiveHistories(EGS_Float effective_histories);
 
    int outputResults(string root_name);
 
    virtual string getInfo() const {
        return getTitle();
    }
 
    bool isValid() const {
        return valid;
    }
 
 
    int outputData(ostream *ofile);
 
    int readData(istream *ifile);
 
    void resetCounter();
 
    int addState(istream &data);
 
 
};
 
class SurfaceCountSpectrum : public BaseSpectrumScorer {
 
 
    void getResult(int bin, EGS_Float &r, EGS_Float &dr);
 
    string getTitle() const {
        return "EGS_Brachy: " + getParticleName() + " counts on surface of source";
    }
 
    string getSubTitle() const {
        return "";
    }
 
    string getYAxisLabel() const {
        return "counts / total counts / MeV\\S-1";
    }
 
    void outputTotal();
 
    string getFileExtension() const {
        return fextension != "" ? fextension : "surfcount";
    }
 
public:
 
 
    SurfaceCountSpectrum(EGS_Input *input, EGS_BaseSource *src, GeomInfo *ginfo, Publisher *publisher):
        BaseSpectrumScorer(input, src, ginfo, publisher) {
 
        publisher->subscribe(this, PARTICLE_ESCAPING_SOURCE);
    };
 
    virtual void score(EB_Message message, void *data=0);
 
 
};
 
 
class EnergyWeightedSurfaceSpectrum : public BaseSpectrumScorer {
 
    void getResult(int bin, EGS_Float &r, EGS_Float &dr);
 
    string getTitle() const {
        return "EGS_Brachy: energy weighted surface " + getParticleName() + " spectrum";
    }
 
    string getSubTitle() const {
        return "";
    }
 
    string getYAxisLabel() const {
        return "(dR/dE)/R / MeV\\S-1";
    }
 
    void outputTotal();
 
    string getFileExtension() const {
        return fextension != "" ? fextension : "ewsurf";
    }
 
public:
 
    EnergyWeightedSurfaceSpectrum(EGS_Input *input, EGS_BaseSource *src, GeomInfo *ginfo, Publisher *publisher):
        BaseSpectrumScorer(input, src, ginfo, publisher) {
 
        publisher->subscribe(this, PARTICLE_ESCAPING_SOURCE);
    };
 
    virtual void score(EB_Message message, void *data=0);
 
 
};
 
class FluenceSpectrumInVoxel : public BaseSpectrumScorer {
 
    EGS_BaseGeometry *geometry;
    EGS_Float region_volume;
    int local_scoring_region;
    int scoring_region;
 
 
    void getResult(int bin, EGS_Float &r, EGS_Float &dr);
 
    string getTitle() const {
        stringstream ss;
        ss << "EGS_Brachy: " << getParticleName() << " fluence spectrum in region ";
        ss << local_scoring_region << " of '"<<geometry->getName()<<"'";
        return ss.str();
    }
 
    string getYAxisLabel() const {
        return "fluence / MeV / cm\\S2";
    }
 
    void outputTotal();
 
    string getFileExtension() const {
        return fextension != "" ? fextension : "voxelflu";
    }
 
public:
 
 
    FluenceSpectrumInVoxel(EGS_Input *input, EGS_BaseSource *src, GeomInfo *ginfo, Publisher *publisher):
        BaseSpectrumScorer(input, src, ginfo, publisher) {
 
        if (particle_type != 0) {
            egsFatal(
                "Energy fluence spectra may only be scored for photons. See:\n\t"
                "https://github.com/nrc-cnrc/EGSnrc/issues/109/\n"
                "for details."
            );
        }
 
        string gname;
        int err = input->getInput("geometry",gname);
        if (err) {
            egsWarning("FluenceSpectrumInVoxel: missing or invalid `geometry` input\n");
            valid =  false;
            return;
        }
 
        geometry = EGS_BaseGeometry::getGeometry(gname);
        if (!geometry) {
            egsWarning("FluenceSpectrumInVoxel: unable to find geometry named `%s` \n", gname.c_str());
            valid =  false;
            return;
        }
 
        err = input->getInput("scoring region", local_scoring_region);
        if (err) {
            egsWarning("FluenceSpectrumInVoxel: no `scoring region` defined. Assuming region 0\n");
            local_scoring_region = 0;
        }
 
        if (scoring_region < 0 || local_scoring_region > geometry->regions()) {
            egsWarning("FluenceSpectrumInVoxel: invalid region `scoring region` defined. Assuming region 0\n");
            local_scoring_region = 0;
        }
 
        GeomRegT geomreg(geometry, local_scoring_region);
        scoring_region = ginfo->localToGlobal(geomreg);
 
        region_volume = geometry->getVolume(local_scoring_region);
 
        publisher->subscribe(this, PARTICLE_TAKING_STEP);
 
    };
 
 
    virtual void score(EB_Message message, void *data=0);
 
 
};
 
 
 
class EnergyFluenceSpectrumInVoxel : public BaseSpectrumScorer {
 
    EGS_BaseGeometry *geometry;
    EGS_Float region_volume;
    int local_scoring_region;
    int scoring_region;
 
 
    void getResult(int bin, EGS_Float &r, EGS_Float &dr);
 
    string getTitle() const {
        stringstream ss;
        ss << "EGS_Brachy: " << getParticleName() << " energy fluence spectrum in region ";
        ss << local_scoring_region << " of '"<<geometry->getName()<<"'";
        return ss.str();
    }
 
    string getYAxisLabel() const {
        return "fluence / MeV / cm\\S2";
    }
 
    void outputTotal();
 
    string getFileExtension() const {
        return fextension != "" ? fextension : "voxelflu";
    }
 
public:
 
 
    EnergyFluenceSpectrumInVoxel(EGS_Input *input, EGS_BaseSource *src, GeomInfo *ginfo, Publisher *publisher):
        BaseSpectrumScorer(input, src, ginfo, publisher) {
 
        if (particle_type != 0) {
            egsFatal(
                "Energy fluence spectra may only be scored for photons. See:\n\t"
                "https://github.com/nrc-cnrc/EGSnrc/issues/109/\n"
                "for details."
            );
        }
 
        string gname;
        int err = input->getInput("geometry",gname);
        if (err) {
            egsWarning("EnergyFluenceSpectrumInVoxel: missing or invalid `geometry` input\n");
            valid =  false;
            return;
        }
 
        geometry = EGS_BaseGeometry::getGeometry(gname);
        if (!geometry) {
            egsWarning("EnergyFluenceSpectrumInVoxel: unable to find geometry named `%s` \n", gname.c_str());
            valid =  false;
            return;
        }
 
        err = input->getInput("scoring region", local_scoring_region);
        if (err) {
            egsWarning("EnergyFluenceSpectrumInVoxel: no `scoring region` defined. Assuming region 0\n");
            local_scoring_region = 0;
        }
 
        if (scoring_region < 0 || local_scoring_region > geometry->regions()) {
            egsWarning("EnergyFluenceSpectrumInVoxel: invalid region `scoring region` defined. Assuming region 0\n");
            local_scoring_region = 0;
        }
 
        GeomRegT geomreg(geometry, local_scoring_region);
        scoring_region = ginfo->localToGlobal(geomreg);
 
        region_volume = geometry->getVolume(local_scoring_region);
 
        publisher->subscribe(this, PARTICLE_TAKING_STEP);
 
    };
 
 
    virtual void score(EB_Message message, void *data=0);
 
 
};