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

Go to the documentation of this file.

/*
################################################################################
#
#  egs_brachy egs_brachy.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
#
################################################################################
*/
 
#ifndef EGS_BRACHY_
#define EGS_BRACHY_
 
// stdlib dependencies
#include <map>
#include <set>
#include <stack>
#include <cstdlib>
 
// egs++ dependencies
#include "egs_advanced_application.h"
#include "egs_functions.h"
#include "egs_input.h"
#include "egs_interface2.h"
#include "egs_interpolator.h"
#include "egs_alias_table.h"
#include "egs_rndm.h"
#include "egs_run_control.h"
#include "egs_scoring.h"
#include "egs_transformations.h"
#include "egs_autoenvelope/egs_autoenvelope.h" // required for Superposition mode
#include "egs_rz/egs_rz.h"
#include "egs_spheres/egs_spheres.h"
 
// local dependencies
#include "pubsub.h"
#include "ginfo.h"
#include "muen.h"
#include "phantom.h"
#include "eb_volcor.h"
#include "spec_scoring.h"
#include "recycle.h"
#include "latch.h"
#include "phsp.h"
#include "timing.h"
 
 
#define PRINT_PARTICLE(P) cout << "\nParticle in reg "<<P.ir << " at "<<P.x.x << " "<<P.x.y<<" "<<P.x.z << " wt "<<P.wt<< " E "<<P.E << " q "<<P.q << " latch "<<P.latch<<"\n";
#define PRINT_PARTICLE_WITH_DIR(P) cout << "\nParticle in reg "<<P.ir << " at "<<P.x.x << " "<<P.x.y<<" "<<P.x.z << " dir "<<P.u.x << " "<<P.u.y<<" "<<P.u.z << " wt "<<P.wt<< " E "<<P.E << " q "<<P.q << " latch "<<P.latch<<"\n";
#define NUM_STUCK_STEPS 1000
#define SAME_POSITION_TOLERANCE 1E-10
#define EB_EPSILON 1E-10
 
class EB_UniformRunControl;
 
class APP_EXPORT EB_Application : public EGS_AdvancedApplication {
 
    enum RunMode {
        RM_NORMAL, 
        RM_SUPERPOSITION, 
        RM_VC_ONLY  
    };
 
 
    RunMode run_mode; 
    string run_mode_name;
 
    EnergyScoringStats *escoring; 
    vector<BaseSpectrumScorer *> spectrum_scorers;
    //SpectrumScoringControl *spec_scoring;
 
    RecycleOpts *recycle_opts;
    bool single_generator;
    bool is_phsp_source;
 
    vector<EGS_Float> source_weights;
    int active_source;
 
    PHSPControl *phsp;
 
    bool score_tlen;   
    bool score_edep;   
    bool score_scat;   
 
    string output_egsdat_format; 
 
    bool output_3ddose_files; 
    string output_dose_format; 
    bool output_egsphant; 
    string output_egsphant_format; 
 
    bool output_voxinfo; 
    string output_voxinfo_format; 
 
    vector<string> output_volcor_phantoms; 
    string output_volcor_format; 
 
    int record_n_init; 
    vector<EGS_Vector> p_init_locs;
 
    EGS_Vector last_position;
    EGS_Float last_R;
    EGS_Float cur_R;
    int steps_at_same_loc;
    EGS_I64 n_stuck;
 
    EGS_BaseGeometry *source_envelope_geom;  
    EGS_ASwitchedEnvelope *superpos_geom;  
 
    vector<EB_Phantom *> phantom_geoms; 
 
    vector<EGS_AffineTransform *> source_transforms; 
    EGS_AffineTransform *base_transform; 
    EGS_AffineTransform *base_transform_inv; 
 
    map<int, EGS_Interpolator *> media_muen; 
    map<string, string> media_muen_names;
 
 
    bool do_brem_split;
    int nbr_split; 
 
    bool do_bcse;
    int bcse_med_num;
    EGS_Float bcse_factor;
    static const EGS_Float DEFAULT_BCSE_FACTOR;
 
    EGS_Float flu_cutoff; 
    EGS_Float source_ecut;  
    EGS_Float source_pcut;  
    EGS_Float global_ecut;  
    EGS_Float global_pcut;  
 
    bool global_i_do_rr;  
    EGS_Float global_e_max_rr;  
 
    bool source_i_do_rr;  
    EGS_Float source_e_max_rr;  
 
    GeomInfo ginfo; 
 
    EB_Phantom *gcr_phantom;  
    int gcr_phantom_reg; 
 
    map<string, vector<int> > extra_scoring_reg;
    map<string, vector<EGS_Float> > extra_scoring_vols;
    map<string, vector<EGS_Float> > extra_scoring_mass;
    map<string, EGS_ScoringArray* > extra_scoring_doses;
    map<string, EGS_ScoringArray* > extra_scoring_doses_edep;
 
    ebvolcor::Results  source_vc_results;  
    ebvolcor::Results  gen_vc_results; 
    ebvolcor::FileResults  file_vc_results; 
 
    Publisher pevent_pub; 
 
    EB_TimingTree timing_blocks; 
 
    map<int, EGS_I64> steps_in_sources;
    map<int, EGS_I64> steps_in_phantoms;
    map<int, EGS_I64> steps_in_other;
 
    ogzstream *gz_data_out;  
    igzstream *gz_data_in;  
 
    static string revision;    
 
 
    int createPhantoms();
 
 
    int correctVolumes();
 
    /* \brief check sources for overlap using a simple MC routine */
    int checkSourceOverlaps(EGS_Input *);
 
    void initTrackLengthScoring(EGS_Input *);
 
    void initMuenData(EGS_Input *);
 
    void initOutputFiles(EGS_Input *);
 
    void initPHSPScoring(EGS_Input *);
 
    void initEDepScoring(EGS_Input *);
 
    void initScatScoring(EGS_Input *);
 
    void clearAusgabCalls();
 
    void enableAusgabCalls(int ncalls, AusgabCall calls[]);
 
 
    void initGCRScoring(EGS_Input *);
 
    void initAusgabCalls();
 
    void initDoseScaling(EGS_Input *);
 
    void initXCCScaling(EGS_Input *);
 
    void initSpectrumScoring(EGS_Input *);
 
    int initCrossSections();
 
    void discardTopParticle(int idisc=1);
 
    void calcEffectiveHistories();
 
    bool isStuck();
 
public:
 
 
    EB_Application(int argc, char **argv) :
        EGS_AdvancedApplication(argc,argv),
        run_mode(RM_NORMAL),
        escoring(),
        recycle_opts(0),
        active_source(0),
        phsp(0),
        score_edep(0),
        score_scat(0),
        output_3ddose_files(true),
        record_n_init(0),
        n_stuck(0),
        superpos_geom(0),
        phantom_geoms(),
        do_brem_split(false),
        nbr_split(1),
        do_bcse(false),
        bcse_factor(1),
        flu_cutoff(0),
        ginfo(),
        gcr_phantom(0), gcr_phantom_reg(0),
        gz_data_out(0),
        gz_data_in(0),
        nsources(0),
        effective_histories(0) {
 
        for (int i=-1; i < 2; i++) {
            steps_in_sources[i] = 0;
            steps_in_phantoms[i] = 0;
            steps_in_other[i] = 0;
        }
 
        pevent_pub.subscribe(&latch_control, PARTICLE_ESCAPING_SOURCE);
        pevent_pub.subscribe(&latch_control, PARTICLE_ESCAPED_SOURCE);
    }
 
    ~EB_Application() {
 
        for (vector<EB_Phantom>::size_type pp=0; pp < phantom_geoms.size(); pp++) {
            delete phantom_geoms[pp];
        }
        phantom_geoms.clear();
 
        map<int, EGS_Interpolator *>::iterator mit = media_muen.begin();
        for (; mit != media_muen.end(); mit++) {
            delete mit->second;
        }
        media_muen.clear();
 
 
        for (vector<EGS_AffineTransform *>::size_type st=0; st < source_transforms.size(); st++) {
            delete source_transforms[st];
        }
        source_transforms.clear();
 
        for (map<string, EGS_ScoringArray* >::iterator sa=extra_scoring_doses.begin(); sa != extra_scoring_doses.end(); sa++) {
            delete sa->second;
        }
        extra_scoring_doses.clear();
 
        for (map<string, EGS_ScoringArray* >::iterator sa=extra_scoring_doses_edep.begin(); sa != extra_scoring_doses_edep.end(); sa++) {
            delete sa->second;
        }
        extra_scoring_doses_edep.clear();
 
        if (phsp) {
            delete phsp;
        }
 
        if (escoring) {
            delete escoring;
        }
 
        vector<BaseSpectrumScorer *>::iterator it = spectrum_scorers.begin();
        for (; it != spectrum_scorers.end(); it++) {
            delete *it;
        }
        spectrum_scorers.clear();
 
        if (recycle_opts) {
            delete recycle_opts;
        }
 
        if (base_transform_inv) {
            delete base_transform_inv;
        }
 
        if (gz_data_in) {
            delete gz_data_in;
        }
 
        if (gz_data_out) {
            delete gz_data_out;
        }
 
    };
 
    int nsources; 
 
    EGS_Float effective_histories;
 
    Latch latch_control;
 
    void describeUserCode() const;
 
    void describeSimulation();
 
    void printIncludedFiles();
 
    int initSimulation();
 
    int initRunControl();
 
    int initRunMode();
 
    int initScoring();
 
    int initVarianceReduction();
 
    int initRussianRoulette(EGS_Input *);
 
    int initBCSE(EGS_Input *);
 
    int ausgab(int iarg);
 
    int outputData();
 
    int outputDataHelper(ostream *);
 
    int egsApplicationOutputData(ostream *);
 
    int egsAdvApplicationOutputData(ostream *);
 
    int egsBrachyOutputData(ostream *);
 
    int readData();
 
    int readDataHelper(istream *);
 
    int egsApplicationReadData(istream *);
 
    int egsAdvApplicationReadData(istream *);
 
    int egsBrachyReadData(istream *);
 
    int combineResults();
 
    void resetCounter();
 
    int addState(istream &data);
 
 
    /* \brief output any results to console */
    void outputResults();
 
 
    void getCurrentResult(double &sum, double &sum2, double &norm,
                          double &count);
 
    EB_Phantom* getPhantomByName(string name);
 
    virtual void startNewParticle();
    virtual void enterNewRegion();
 
    virtual int runSimulation();
 
    string getOutputVolcorFormat() {
        return output_volcor_format;
    }
 
    int howManyJobsDone();
 
protected:
 
    void addRecycledParticlesToStack(EGS_Particle *p, bool new_hist=false);
    void copyParticleToSourceLoc(EGS_Particle *p, int source, bool kill_orig, bool rotate, EGS_Float new_wt);
 
    void doPhotonSplitting(int);
 
 
    /* \brief initialize a single shower history.
     *
     * Particle is initiated from source and then transformed to location
     * of next source */
    int simulateSingleShower();
    int startNewShower();
 
    int initGeometry();
 
 
    /* \brief override default initSource to allow us to get source position weighting */
    int initSource();
 
 
    int initSourceTransforms();
 
    /* \brief take an input containging multiple `transformation` input blocks
     * and return a vector of pointers to the generated transforms.
     */
    vector<EGS_AffineTransform *> createTransforms(EGS_Input *input);
 
};
 
class EB_UniformRunControl : public EGS_RunControl {
 
public:
 
    EB_UniformRunControl(EB_Application *app, string egsdat_format);
    ~EB_UniformRunControl() {};
 
    void  describeRCO();
 
    int  startSimulation();
 
    int  finishSimulation();
 
protected:
 
    EB_Application *app; 
 
    string output_egsdat_format; 
    
    int milliseconds;   // time interval for checking
    // if all jobs finished (default 1000 ms)
 
    int check_intervals;// Number of intervals to check
    // if all jobs done (default 5)
 
    int    njob;
    int    npar;
    int    ipar;
    int    ifirst;
    bool   check_egsdat;// If true, and a 'watcher' job, produce intermediate results
    bool   watcher_job; // If true, job is a 'watcher'
};
 
#endif