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 0998343e8..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 @@ -68,22 +68,22 @@ void EGS_RadiativeSplitting::setApplication(EGS_Application *App) { char buf[32]; // Set EGSnrc internal UBS + RR - if ( i_play_RR ){ - app->setRussianRoulette(nsplit); - i_play_RR = true; + 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; + 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 ( i_play_RR ){ + if (i_play_RR) { description +="\n - Splitting radiative events in "; sprintf(buf,"%d",nsplit); description += buf; 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 58dd84a84..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 @@ -109,26 +109,26 @@ class EGS_RADIATIVE_SPLITTING_EXPORT EGS_RadiativeSplitting : public EGS_AusgabO /*! 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; + 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; }; 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 ){ + 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 { @@ -143,27 +143,27 @@ class EGS_RADIATIVE_SPLITTING_EXPORT EGS_RadiativeSplitting : public EGS_AusgabO 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); + 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); + /* 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); + /* 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); } diff --git a/HEN_HOUSE/egs++/egs_advanced_application.cpp b/HEN_HOUSE/egs++/egs_advanced_application.cpp index de699caeb..350e99072 100644 --- a/HEN_HOUSE/egs++/egs_advanced_application.cpp +++ b/HEN_HOUSE/egs++/egs_advanced_application.cpp @@ -1227,67 +1227,82 @@ void EGS_AdvancedApplication::setRadiativeSplitting(const EGS_Float &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; + 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; + 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; +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; + 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; } //************************************************************ diff --git a/HEN_HOUSE/egs++/egs_advanced_application.h b/HEN_HOUSE/egs++/egs_advanced_application.h index 4b560b091..683d0cd71 100644 --- a/HEN_HOUSE/egs++/egs_advanced_application.h +++ b/HEN_HOUSE/egs++/egs_advanced_application.h @@ -236,10 +236,10 @@ class APP_EXPORT EGS_AdvancedApplication : public EGS_Application { //************************************************ /* Turn ON/OFF EGSnrc internal radiative splitting (UBS) */ - void setRadiativeSplitting( const EGS_Float &nsplit ); + 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 ); + 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 d3fd26d85..59fc1d45e 100644 --- a/HEN_HOUSE/egs++/egs_application.h +++ b/HEN_HOUSE/egs++/egs_application.h @@ -1155,7 +1155,7 @@ class EGS_EXPORT EGS_Application { //************************************************ virtual void setRadiativeSplitting(const EGS_Float &nsplit) {}; virtual void setRussianRoulette(const EGS_Float &iSwitchRR) {}; - virtual void splitTopParticleIsotropically(const EGS_Float &fsplit){} + virtual void splitTopParticleIsotropically(const EGS_Float &fsplit) {} //************************************************************ // Utility functions for use with ausgab fluence scoring objects