Fix radiative splitting ao

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<bool> is_sensitive;     // flag scoring regions
-    vector<bool> is_source;        // Flag regions such as brems target or radiactive source
+    vector<bool> is_source;        // Flag regions such as brems target or radioactive source
     // Interacting particles not subjected to classification
     vector <int> f_start, f_stop;  // Markers for group regions input
     vector <int> f_region;         // Input list of scoring regions

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);
 
         //=================================================

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;
 
 };
 

HEN_HOUSE/egs++/egs_advanced_application.cpp

@@ -1213,11 +1213,81 @@ 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
 //************************************************************

HEN_HOUSE/egs++/egs_advanced_application.h

@@ -229,8 +229,16 @@ class APP_EXPORT EGS_AdvancedApplication : public EGS_Application {
 
     /* Needed by some sources */
     EGS_Float getRM();
-    /* Turn ON/OFF radiative splitting */
-    void setRadiativeSplitting(const EGS_Float &nsplit);
+
+    //************************************************
+    // 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:
 

HEN_HOUSE/egs++/egs_application.h

@@ -1148,7 +1148,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