diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_fluence_scoring/egs_fluence_scoring.h b/HEN_HOUSE/egs++/ausgab_objects/egs_fluence_scoring/egs_fluence_scoring.h index 9bf2be808..433b2f9fc 100644 --- a/HEN_HOUSE/egs++/ausgab_objects/egs_fluence_scoring/egs_fluence_scoring.h +++ b/HEN_HOUSE/egs++/ausgab_objects/egs_fluence_scoring/egs_fluence_scoring.h @@ -200,7 +200,7 @@ class EGS_FLUENCE_SCORING_EXPORT EGS_FluenceScoring : public EGS_AusgabObject { /* Regions flags */ vector is_sensitive; // flag scoring regions - vector is_source; // Flag regions such as brems target or radiactive source + vector is_source; // Flag regions such as brems target or radioactive source // Interacting particles not subjected to classification vector f_start, f_stop; // Markers for group regions input vector f_region; // Input list of scoring regions diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/egs_radiative_splitting.cpp b/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/egs_radiative_splitting.cpp index fbcd05b52..613b70c9f 100644 --- a/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/egs_radiative_splitting.cpp +++ b/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/egs_radiative_splitting.cpp @@ -67,15 +67,31 @@ void EGS_RadiativeSplitting::setApplication(EGS_Application *App) { char buf[32]; - // Set EGSnrc internal radiative splitting number . - app->setRadiativeSplitting(nsplit); + // Set EGSnrc internal UBS + RR + if (i_play_RR) { + app->setRussianRoulette(nsplit); + i_play_RR = true; + } + // Set EGSnrc internal radiative splitting number. + else if (nsplit > 0) { + app->setRadiativeSplitting(nsplit); + i_play_RR = false; + } description = "\n===========================================\n"; description += "Radiative splitting Object ("; description += name; description += ")\n"; description += "===========================================\n"; - if (nsplit > 1) { + if (i_play_RR) { + description +="\n - Splitting radiative events in "; + sprintf(buf,"%d",nsplit); + description += buf; + description +="\n - Play RR with higher order e-/e+ with probability 1/"; + sprintf(buf,"%d\n\n",nsplit); + description += buf; + } + else if (nsplit > 1) { description +="\n - Splitting radiative events in "; sprintf(buf,"%d\n\n",nsplit); description += buf; @@ -104,6 +120,7 @@ extern "C" { } EGS_Float nsplit = 1.0; + /*! Switch for splitting + RR. Negative nsplit value switches OFF RR. */ int err = input->getInput("splitting",nsplit); //================================================= diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/egs_radiative_splitting.h b/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/egs_radiative_splitting.h index 7de1b99bf..cff1c2129 100644 --- a/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/egs_radiative_splitting.h +++ b/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/egs_radiative_splitting.h @@ -106,21 +106,77 @@ class EGS_RADIATIVE_SPLITTING_EXPORT EGS_RadiativeSplitting : public EGS_AusgabO ~EGS_RadiativeSplitting(); void setApplication(EGS_Application *App); - + /*! Switch for splitting + RR. Negative nsplit value switches OFF RR. */ void setSplitting(const int &n_s) { nsplit = n_s; + if (nsplit < 0) { + nsplit *= -1; + i_play_RR = false; + } + else if (nsplit > 1) { + i_play_RR = true; + } + /* Avoid zero division. A zero value turns off brems */ + wthin = nsplit ? 1./nsplit : 1.0; }; - int processEvent(EGS_Application::AusgabCall iarg) { - return 0; + bool needsCall(EGS_Application::AusgabCall iarg) const { + if ( + iarg == EGS_Application::BeforeBrems || + iarg == EGS_Application::BeforeAnnihFlight || + iarg == EGS_Application::BeforeAnnihRest || + iarg == EGS_Application::AfterBrems || + iarg == EGS_Application::AfterAnnihFlight || + iarg == EGS_Application::AfterAnnihRest || + iarg == EGS_Application::FluorescentEvent) { + return true; + } + else { + return false; + } }; - int processEvent(EGS_Application::AusgabCall iarg, int ir) { + + int processEvent(EGS_Application::AusgabCall iarg) { + + /* A fat particle's weight is larger than a thin particle's max weight */ + bool is_phat = (app->top_p.wt - wthin) > epsilon; + bool is_primary = app->top_p.latch == 0 ? true : false; + + /* Split radiative events ONLY for primary and fat electrons */ + if (iarg == EGS_Application::BeforeBrems || + iarg == EGS_Application::BeforeAnnihFlight || + iarg == EGS_Application::BeforeAnnihRest && + (is_primary || is_phat)) { + app->setRadiativeSplitting(nsplit); + } + /* Avoids higher order splitting of radiative events */ + else if (iarg == EGS_Application::AfterBrems || + iarg == EGS_Application::AfterAnnihFlight || + iarg == EGS_Application::AfterAnnihRest) { + app->setRadiativeSplitting(1); + app->setLatch(app->getNpOld()+1,1); + } + /* Fluorescent photons created by charged particles surviving RR + when radiative splitting ON should be split to avoid having heavy photons. + This should happen in EGSnrc, but it is not implemented yet, so do it here! + Note that when this is implemented in EGSnrc, the weight check will make sure + photons aren't split again! + */ + if (iarg == EGS_Application::FluorescentEvent && is_phat && nsplit > 1) { + app->splitTopParticleIsotropically(nsplit); + } + + return 0; }; protected: - /* Maximum splitting limited to 2,147,483,647 */ + /* Max weight of thin particles */ + EGS_Float wthin; + /* Maximum splitting limited to 2,147,483,647. Negative value switches OFF RR. */ int nsplit; + /* Switch for Russian Roulette */ + bool i_play_RR; }; diff --git a/HEN_HOUSE/egs++/egs_advanced_application.cpp b/HEN_HOUSE/egs++/egs_advanced_application.cpp index e229d30ea..350e99072 100644 --- a/HEN_HOUSE/egs++/egs_advanced_application.cpp +++ b/HEN_HOUSE/egs++/egs_advanced_application.cpp @@ -1215,11 +1215,96 @@ EGS_Float EGS_AdvancedApplication::getPcut() { EGS_Float EGS_AdvancedApplication::getRM() { return the_useful->rm; } -// Turns ON/OFF radiative splitting + +//************************************************************ +// Utility functions for fluence scoring objects +//************************************************************ + +// Turns ON/OFF EGSnrc internal radiative splitting (UBS) void EGS_AdvancedApplication::setRadiativeSplitting(const EGS_Float &nsplit) { the_egsvr->nbr_split = nsplit; } +// Turns ON/OFF EGSnrc internal Russian Roulette + UBS +void EGS_AdvancedApplication::setRussianRoulette(const EGS_Float &iSwitchRR) { + if (iSwitchRR > 1.0) { + the_egsvr->i_play_RR = 1; + the_egsvr->prob_RR = 1.0/iSwitchRR; + the_egsvr->nbr_split = iSwitchRR; + } + else { + the_egsvr->i_play_RR = 0; + the_egsvr->nbr_split = 1; + } +} + +// Splits top particle into nsplit particles uniformly in 4Pi +void EGS_AdvancedApplication::splitTopParticleIsotropically(const EGS_Float &fsplit) { + // Reset particle pointer + the_stack->npold = the_stack->np; + /* Initialize local variables */ + int np = the_stack->np-1, + the_latch = the_stack->latch[np], + ir = the_stack->ir[np], + iq = the_stack->iq[np]; + the_stack->wt[np] /= fsplit; + double E = the_stack->E[np]; + EGS_Float x = the_stack->x[np], y = the_stack->y[np], z = the_stack->z[np], + wthin = the_stack->wt[np], dnear = the_stack->dnear[np]; + EGS_Float u,v,w; + /* If fsplit is a non-integer, sample between int(fsplit) and int(split)+1 */ + int nsplit = int(fsplit); + EGS_Float dsplit = fsplit - nsplit; + if (dsplit > 0) { // non-integer splitting number + if (rndm->getUniform() < dsplit) { + ++nsplit; + } + } + for (int i=0; i < nsplit; i++) { + np++; + if (np >= MXSTACK) { + egsFatal("\n\n******************************************\n" + "ERROR: In EGS_AdvancedApplication::splitTopParticleIsotropically() :\n" + "max. stack depth MXSTACK=%d < np=%d\n" + "Stack overflow due to splitting!\n" + "******************************************\n" + ,MXSTACK,np); + } + the_stack->x[np] = x; + the_stack->y[np] = y; + the_stack->z[np] = z; + the_stack->iq[np]= iq; + the_stack->dnear[np] = dnear; + the_stack->latch[np] = the_latch; + the_stack->ir[np]= ir; + the_stack->E[np] = E; + the_stack->wt[np]=wthin; + // Particles isotropically distributed in space + w = 2*rndm->getUniform()-1; + EGS_Float sinz = 1-w*w; + if (sinz > epsilon) { + sinz = sqrt(sinz); + EGS_Float cphi, sphi; + EGS_Float phi = 2*M_PI*rndm->getUniform(); + cphi = cos(phi); + sphi = sin(phi); + u = sinz*cphi; + v = sinz*sphi; + } + else { + u = 0; + v = 0; + } + + the_stack->u[np] = u; + the_stack->v[np] = v; + the_stack->w[np] = w; + + } + + the_stack->np = np+1; +} + //************************************************************ // Utility functions for fluence scoring objects //************************************************************ diff --git a/HEN_HOUSE/egs++/egs_advanced_application.h b/HEN_HOUSE/egs++/egs_advanced_application.h index c69373424..683d0cd71 100644 --- a/HEN_HOUSE/egs++/egs_advanced_application.h +++ b/HEN_HOUSE/egs++/egs_advanced_application.h @@ -230,8 +230,16 @@ class APP_EXPORT EGS_AdvancedApplication : public EGS_Application { /* Needed by some sources */ EGS_Float getRM(); - /* Turn ON/OFF radiative splitting */ + + //************************************************ + // For use with ausgab radiative splitting objects + //************************************************ + + /* Turn ON/OFF EGSnrc internal radiative splitting (UBS) */ void setRadiativeSplitting(const EGS_Float &nsplit); + /* Turn ON/OFF EGSnrc internal Russian Roultette + UBS */ + void setRussianRoulette(const EGS_Float &iSwitchRR); + void splitTopParticleIsotropically(const EGS_Float &fsplit); protected: diff --git a/HEN_HOUSE/egs++/egs_application.h b/HEN_HOUSE/egs++/egs_application.h index 960ab4c20..59fc1d45e 100644 --- a/HEN_HOUSE/egs++/egs_application.h +++ b/HEN_HOUSE/egs++/egs_application.h @@ -1149,7 +1149,13 @@ class EGS_EXPORT EGS_Application { virtual EGS_Float getRM() { return -1.0; }; + + //************************************************ + // For use with ausgab radiative splitting objects + //************************************************ virtual void setRadiativeSplitting(const EGS_Float &nsplit) {}; + virtual void setRussianRoulette(const EGS_Float &iSwitchRR) {}; + virtual void splitTopParticleIsotropically(const EGS_Float &fsplit) {} //************************************************************ // Utility functions for use with ausgab fluence scoring objects