diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/Makefile b/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/Makefile index ab04c0b7b..36dfb6ac2 100644 --- a/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/Makefile +++ b/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/Makefile @@ -23,7 +23,7 @@ # # Author: Ernesto Mainegra-Hing, 2018 # -# Contributors: +# Contributors: Blake Walters, 2021 # ############################################################################### @@ -41,8 +41,14 @@ extra_dep = $(addprefix $(DSOLIBS), $(my_deps)) include $(SPEC_DIR)egspp_libs.spec +INC2 = -I$(IEGS2) -I..$(DSEP).. + $(make_depend) +clean: + $(REMOVE) $(ABS_DSO)$(library)* + test: @echo "common_h2: $(common_h2)" + @echo "INC2: $(INC2)" @echo "extra_dep: $(extra_dep)" diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/array_sizes.h b/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/array_sizes.h new file mode 100644 index 000000000..ff8d13668 --- /dev/null +++ b/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/array_sizes.h @@ -0,0 +1,52 @@ +/* +############################################################################### +# +# EGSnrc radiative splitting ausgab object array sizes headers +# Copyright (C) 2015 National Research Council Canada +# +# This file is part of EGSnrc. +# +# EGSnrc 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. +# +# EGSnrc 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. +# +# You should have received a copy of the GNU Affero General Public License +# along with EGSnrc. If not, see . +# +############################################################################### +# +# Author: Iwan Kawrakow, 2005 +# +# Contributors: +# +############################################################################### +# +# Defines he maximum number of media (MXMED) and the maximum number of +# particles on the stack (MXSTACK). This file gets included by the egsnrc +# fortran subroutines (egsnrc_$my_machine.F), the base application +# (egs_simple_application.cpp or egs_advanced_application.cpp in +# $HEN_HOUSE/egs++), and possibly the user code, if it uses the particle +# stack or one of the structures that depends on the maximum number of media. +# +############################################################################### +*/ + + +#ifndef ARRAY_SIZES_ +#define ARRAY_SIZES_ +#define MXMED 10 +#define MXSTACK 1000000 +#define MXBRES 100 +#define MXBRXS 50 +#define MXEL 50 +#define MXPWR2I 50 +#define MXGE 2000 +#define MXRAYFF 100 +#define MXSGE 400 +#endif 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..538d60508 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 @@ -23,7 +23,7 @@ # # Author: Ernesto Mainegra-Hing, 2018 # -# Contributors: +# Contributors: Blake Walters, 2021 # ############################################################################### # @@ -31,7 +31,7 @@ # # TODO: # -# - Add directional radiative splitting (DRS) +# - Testing/debugging directional radiative splitting (DRS) # ############################################################################### */ @@ -45,10 +45,82 @@ #include #include #include +#include #include "egs_radiative_splitting.h" #include "egs_input.h" #include "egs_functions.h" +#include "egs_transformations.h" +#include "egs_math.h" + +//local structures required +struct DBS_Aux { + int np; + EGS_Float E; + EGS_Float Z13; + EGS_Float le; +}; + +//local functions +EGS_Float dbs_func_KM(EGS_Float x, void *data) { + DBS_Aux *the_aux = (DBS_Aux *)data; + EGS_Float Z13 = the_aux->Z13; + int np = the_aux->np; + EGS_Float E = the_aux->E; + EGS_Float k = (E-1)*x; + EGS_Float Ep = E - k; + EGS_Float r = Ep/E; + EGS_Float delta = k/(2*E*Ep); + delta *= delta; + EGS_Float r2p1 = 1+r*r, rp12 = r2p1 + 2*r; + EGS_Float beta = sqrt((E-1)*(E+1))/E; + EGS_Float y2min = 0, y2max = 2*beta*(1+beta)*E*E; + EGS_Float y2maxi = 1/y2max; + EGS_Float a = sqrt(1+y2max); + EGS_Float deta = 1./np; + double sum = 0; + for(int j=0; jZ13; + int np = the_aux->np; + EGS_Float E = the_aux->E; + EGS_Float k = E - exp(u*the_aux->le); + EGS_Float Ep = E - k; + EGS_Float r = Ep/E; + EGS_Float delta = k/(2*E*Ep); + delta *= delta; + EGS_Float r2p1 = 1+r*r, rp12 = r2p1 + 2*r; + EGS_Float beta = sqrt((E-1)*(E+1))/E; + EGS_Float y2min = 0, y2max = 2*beta*(1+beta)*E*E; + EGS_Float y2maxi = 1/y2max; + EGS_Float a = sqrt(1+y2max); + EGS_Float deta = 1./np; + double sum = 0; + for(int j=0; jsetRadiativeSplitting(nsplit); + if (( split_type == DRS || split_type == DRSf ) && app->getIbrdst()) + { + initSmartKM(app->getEmax()*1.05); + } + description = "\n===========================================\n"; description += "Radiative splitting Object ("; description += name; @@ -77,18 +154,1510 @@ void EGS_RadiativeSplitting::setApplication(EGS_Application *App) { description += "===========================================\n"; if (nsplit > 1) { description +="\n - Splitting radiative events in "; - sprintf(buf,"%d\n\n",nsplit); + sprintf(buf,"%d",nsplit); description += buf; + if ( split_type == URS ) { + description +=" using URS\n\n"; + } + else if ( split_type == DRS || split_type == DRSf) { + description +=" using DRS with the following parameters:\n"; + description +=" fs = "; + sprintf(buf,"%g",fs); + description += buf; + description += " cm\n"; + description +=" ssd = "; + sprintf(buf,"%g",ssd); + description += buf; + description += " cm\n\n"; + if (ireg_esplit.size()>0) { + //we are splitting e-/e+ + description += " phat e-/e+ will be split "; + sprintf(buf,"%d",nsplit); + description += buf; + description += " times on entering the following regions: "; + for (int i=0; igetRotation().xx(),T->getRotation().xy(),T->getRotation().xz()); + description += buf; + sprintf(buf," %.5f %.5f %.5f\n",T->getRotation().yx(),T->getRotation().yy(),T->getRotation().yz()); + description += buf; + sprintf(buf," %.5f %.5f %.5f\n",T->getRotation().zx(),T->getRotation().zy(),T->getRotation().zz()); + description += buf; + description += " translation vector = "; + sprintf(buf,"%.5f %.5f %.5f\n\n",T->getTranslation().x,T->getTranslation().y,T->getTranslation().z); + description += buf; + } + } + else { + description +=" using an unknown splitting algorithm? :-("; + } } - else if (nsplit == 1) { + if (nsplit == 1) { description +="\n - NO radiative splitting"; } else { - description +="\n - BEWARE: Turning OFF radiative events !!!"; + description +="\n - BEWARE: Turning OFF EGSnrc radiative events !!!"; } + description += "\n===========================================\n\n"; } +void EGS_RadiativeSplitting::initDBS(const EGS_Float field_rad, const EGS_Float field_ssd, const vector splitreg, const int irad, const EGS_Float zrr, const EGS_AffineTransform *t) { + fs = field_rad; + ssd = field_ssd; + ireg_esplit = splitreg; + irad_esplit = irad; + zrr_esplit = zrr; + if (T) { + //delete T; + T = 0; + } + if (t) { + T = new EGS_AffineTransform(*t); + if (!T) { + egsWarning("Invalid transform for DBS radiative splitting cone."); + } + } +} + +//at least try to get brems working with this +int EGS_RadiativeSplitting::doInteractions(int iarg, int &killed) +{ + + //data for initiating particle + int np = app->getNp(); + EGS_Particle pi = app->getParticleFromStack(np); + //DBS calculations assume splitting field is centred on Z-axis and ssd is defined relative to Z=0 + inverseTransformP(pi,T); + EGS_Float dneari = app->getDnear(np); + int latch = pi.latch; + + int check = 0; //set to 1 to return to shower + + killed = 0; + + if( iarg > EGS_Application::AfterTransport && pi.x.z > ssd ) { + //particle is past ssd, no splitting + app->setRadiativeSplitting(1); + return 0; + } + + //use bit 0 to mark phat particles + //seems like a temporary solution + int is_fat = (latch & (1 << 0)); + + //below was used for debugging but could be true in general + //if(pi.wt < 1 && is_fat) exit(1); + + if( iarg == EGS_Application::BeforeTransport) + { + return 0; + } + + if( iarg == EGS_Application::AfterTransport) + { + if (pi.q && pi.u.z > 0 && is_fat && ireg_esplit.size()>0) + { + //see if charged particle has entered a splitting region + auto it = std::find(ireg_esplit.begin(),ireg_esplit.end(),pi.ir); + if (it != ireg_esplit.end()) + { + //split charged particle nsplit times and radially redistributed (if required) + EGS_Float wt = pi.wt/nsplit; + latch = latch & ~(1 << 0); + EGS_Float ang_dbs = 0; + if (irad_esplit) + { + ang_dbs = 2*M_PI/nsplit; + } + //delete the charged particle and replace with nsplit particles (possibly redistributed + //about the Z-axis in untransformed space) + app->deleteParticleFromStack(np); + for (int i=0; iaddParticleToStack(p,dneari); + } + } + } + return 0; + } + + if( iarg == EGS_Application::BeforeBrems ) { + double E = pi.E; + EGS_Float wt = pi.wt; + if( is_fat ) { + //clear bit 0 + latch = latch & ~(1 << 0); + //reset latch value of top particle + app->setLatch(latch); + //is the next line necessary? + app->setRadiativeSplitting(nsplit); + int npold = np; + int res = doSmartBrems(); + if( res ) { + app->callBrems(); + int nstart = np+1, aux=0; + killThePhotons(fs,ssd,nsplit,nstart,aux); + } + //we need to relable the interacting e- as fat + //at this point np holds the e- + EGS_Particle p = app->getParticleFromStack(npold); + latch = latch | (1 << 0); + p.latch = latch; + app->updateParticleOnStack(npold,p,dneari); + } + else { + app->setRadiativeSplitting(1); + app->callBrems(); + int nstart = np+1, aux=0; + killThePhotons(fs,ssd,nsplit,nstart,aux); + } + check = 1; + } + else if( iarg == EGS_Application::BeforeAnnihFlight ) { + if( is_fat ) { + //figure out what to do with the extra stack + //the_extra_stack->iweight[np]=1; + //set interacting particle to nonphat so this is + //passed on to resultant photons + //TODO: figure out a better way to do this that does not use latch + latch = latch & ~(1 << 0); + app->setLatch(latch); + app->setRadiativeSplitting(nsplit); + } + else app->setRadiativeSplitting(1); + app->callAnnih(); + int nstart=np+1,aux=0; + killThePhotons(fs,ssd,nsplit,nstart,aux); + check = 1; + } + else if( iarg == EGS_Application::BeforeAnnihRest ) { + if( is_fat ) { + //app->setNpold(np);not sure why this is necessary--unless it is intended to ensure + //overwriting of the annihilated positron, in which case, we should + //actually be setting npold=np-1 + app->setRadiativeSplitting(nsplit);//this also seems like a formality + int nsamp = 2*nsplit; + //label photons as nonphat + latch = latch & ~(1 << 0); + app->setLatch(latch); + uniformPhotons(nsamp,nsplit,fs,ssd,app->getRM()); + } + else { + app->setRadiativeSplitting(1); + app->callAnnihAtRest(); + int nstart=np+1,aux=0; + killThePhotons(fs,ssd,nsplit,nstart,aux); + } + check = 1; + } + if( iarg == EGS_Application::BeforeCompton || + iarg == EGS_Application::BeforePair || + iarg == EGS_Application::BeforePhoto || + iarg == EGS_Application::BeforeRayleigh ) + { + + //if e- splitting is off or Z <= the Russian Roulette plane, only allow phat photons to undergo interaction + //unless it is Rayleigh or bound compton + if(!is_fat && (iarg == EGS_Application::BeforePhoto || iarg == EGS_Application::BeforePair || + (iarg == EGS_Application::BeforeCompton && !app->getIbcmp())) && + (ireg_esplit.size()==0 || pi.x.z <= zrr_esplit)) + { + if (app->getRngUniform()*nsplit > 1) + { + //send interacting particle back with zero wt + pi.wt = 0; + app->updateParticleOnStack(np,pi,dneari); + return 0; + } + else + { + //keep the photon, increase weight and label as phat + pi.wt = pi.wt*nsplit; + pi.latch = pi.latch | (1 << 0); + //copy pi because we are transforming it before putting it back on the stack + EGS_Particle p = pi; + transformP(p,T); + app->updateParticleOnStack(np,p,dneari); + is_fat = 1; + } + } + + int nint = is_fat ? nsplit : 1; + + if(iarg == EGS_Application::BeforePair) + { + //only split phat interactions if e-/e+ splitting is on and Z > zrr_esplit + if (ireg_esplit.size() > 0 && pi.x.z > zrr_esplit) + { + latch = latch & ~(1 << 0); + EGS_Particle p = pi; + p.latch = latch; + p.wt = p.wt/nint; + //Do transform here. Seems more efficient than inside the loop below. + transformP(p,T); + app->deleteParticleFromStack(np); + for (int i=0; iaddParticleToStack(p,dneari); + app->callPair(); + } + } + else + { + //let EGSnrc take care of the interaction + app->callPair(); + } + } + if(iarg == EGS_Application::BeforePhoto) + { + //only split phat interactions if e-/e+ splitting is on and Z > zrr_esplit + if (ireg_esplit.size() > 0 && pi.x.z > zrr_esplit) + { + //label interacting photon as non-phat and reduce weight (if phat) + latch = latch & ~(1 << 0); + EGS_Particle p = pi; + p.latch = latch; + p.wt = p.wt/nint; + transformP(p,T); + app->deleteParticleFromStack(np); + for (int i=0; iaddParticleToStack(p,dneari); + app->callPhoto(); + EGS_Particle pj = app->getParticleFromStack(app->getNp()); + } + } + else + { + app->callPhoto(); + } + } + if(iarg == EGS_Application::BeforeCompton) + { + int npold = np; + if(is_fat && !app->getIbcmp() && (ireg_esplit.size()==0 || pi.x.z <= zrr_esplit) ) + { + //label as nonphat to be passed on to descendents + latch = latch & ~(1 << 0); + EGS_Particle p = pi; + p.latch = latch; + transformP(p,T); + app->updateParticleOnStack(np,p,dneari); + doSmartCompton(nint); + } + else //straight-up compton + { + //label interacting photon as non-phat and reduce weight (if phat) + latch = latch & ~(1 << 0); + EGS_Particle p = pi; + p.latch = latch; + p.wt = p.wt/nint; + transformP(p,T); + app->deleteParticleFromStack(np); + for (int i=0; iaddParticleToStack(p,dneari); + int nstart = app->getNp(); + app->callCompt(); + int aux = 1; + if (nint == 1 || (ireg_esplit.size()>0 && pi.x.z > zrr_esplit)) + { + //not sure why we do not kill the electron if nint=1, but this is copied from BEAMnrc + aux = 0; + } + killThePhotons(fs,ssd,nsplit,nstart,aux); + } + } + } + else if (iarg == EGS_Application::BeforeRayleigh) + { + //TODO: Put all of this in a doRayleigh function + if (is_fat) { + latch = latch & ~(1 << 0); + } + EGS_Float gle = app->getGle(); + int imed = app->getMedium(pi.ir); + int lgle = app->getLgle(gle,imed); + EGS_Float costhe, sinthe; + //delete the interacting particle because we are going to replace it + app->deleteParticleFromStack(np); + for (int i=0; icallEgsRayleighSampling(imed,p.E,gle,lgle,costhe,sinthe); + //adjust scatter angles and apply to particle + doUphi21(sinthe,costhe,p.u); + //add the particle to the stack + transformP(p,T); + app->addParticleToStack(p,dneari); + //now potentially kill it -- seems like we should kill the particle before adding to the stack + int nstart = app->getNp(), aux=0; + killThePhotons(fs,ssd,nsplit,nstart,aux); + } + } + check = 1; + } + else if(iarg == EGS_Application::FluorescentEvent ) + { + if( is_fat ) { + EGS_Float ener = pi.E; + //label photons as nonphat + latch = latch & ~(1 << 0); + app->setLatch(latch); + uniformPhotons(nsplit,nsplit,fs,ssd,ener); + } + else { + int nstart=np, aux=0; + killThePhotons(fs,ssd,nsplit,nstart,aux); + } + // Do not need to return to shower + check = 2; + } + + /* + //summary debug statement + for(int i=np; i<=app->getNp(); i++) + { + EGS_Particle p = app->getParticleFromStack(i); + int is_fat_test = (p.latch & (1 << 0)); + if (p.x.z < 0) + { + egsInformation("after everything: iarg=%d i=%d np=%d iq=%d E=%g wt=%g is_fat=%d x=%g y=%g z=%g u=%g v=%g w=%g\n",iarg,i,np,p.q,p.E,p.wt,is_fat_test,p.x.x,p.x.y,p.x.z,p.u.x,p.u.y,p.u.z); + } + } + */ + + if( check ) { + if( check == 1 ) { + if( app->getNp() >= app->getMxstack() ) egsFatal("Stack size exceeded " + "in EGS_RadiativeSplitting::doInteractions()\n"); + //set weight of particle on top of stack to 0 + EGS_Particle p = app->getParticleFromStack(app->getNp()); + EGS_Float dnear = app->getDnear(app->getNp()); + p.wt = 0.; + app->addParticleToStack(p,dnear); + } + return 0; + } + + return 1; +} + +int EGS_RadiativeSplitting::doSmartBrems() { + int np = app->getNp(); + int iwt = nsplit; + int nbrspl = iwt; + + //clear stack of brems particles + particle_stack.clear(); + + EGS_Particle pi = app->getParticleFromStack(np); + inverseTransformP(pi,T); + EGS_Float dneari = app->getDnear(np); + + double E = pi.E; + EGS_Float ener = E - app->getRM(); + EGS_Float tau = ener/app->getRM(); + EGS_Float beta = sqrt(tau*(tau+2))/(tau+1); + EGS_Float beta2 = beta*beta; + EGS_Vector x = pi.x; + EGS_Vector u = pi.u; + int irl = pi.ir; + int latch = pi.latch; + + EGS_Float ct_min,ct_max,ro; + getCostMinMax(x,u,ro,ct_min,ct_max); + imed = app->getMedium(irl); + EGS_Float f_max_KM = 1, q_KM, p_KM; + int j_KM; + + //for now do not use doSmartBrems when IBRDST=1 (full KM expression used for brems angular sampling) + //this may be implemented in the future with another DBS option in which we use the implementation from beampp + //For this reason we have not deleted coding for ibrdbst=1 below. + if(app->getIbrdst() == 1) { + return 1; + } + + if(app->getIbrdst() == 1) { + q_KM = a_KM[imed]*log(ener) + b_KM[imed]; + j_KM = (int) q_KM; + q_KM -= j_KM; + p_KM = 1 - q_KM; + f_max_KM = f_KM_max[imed].interpolateFast(j_KM,log(ener)); + } + + EGS_Float w1, cmin, cmax; + if( app->getIbrdst() == 1 ) { + cmin = sqrt(1-beta*ct_min); + cmax = sqrt(1-beta*ct_max); + w1 = (cmin - cmax)*sqrt(1+beta)/(cmax*cmin*((1+beta)*(1+tau)-1)); + } + else w1 = (ct_max - ct_min)/((1-beta*ct_max)*(1-beta*ct_min)*2* + (tau+1)*(tau+1)); + w1 *= f_max_KM; + EGS_Float d = ssd - x.z; + EGS_Float dmin = ro <= fs ? d : sqrt(d*d + (ro-fs)*(ro-fs)); + EGS_Float w2 = fs*fs*d/(2*dmin*dmin*dmin); + bool will_rotate = use_cyl_sym && x.z < zcyls; + + //if( will_rotate ) w2 *= rsamp->getAeff(); + EGS_Float aux; + if( app->getIbrdst() == 1 ) { + w2 *= beta*sqrt(1+beta)/(2*(1-beta*ct_max)*sqrt(1-beta*ct_max)* + ((1+beta)*(tau+1)-1)); + } + else { + aux = (tau+1)*(1-beta*ct_max); + w2 /= (2*aux*aux); + } + w2 *= f_max_KM; + EGS_Float wprob = min(w1,w2); + if( wprob >= 1 && app->getIbrdst() == 1 ) return 1; + int nsample; + if( wprob > 1 ) { + ct_min = -1; + ct_max = 1; + wprob = 1; + } + EGS_Float asample = wprob*nbrspl; + nsample = (int) asample; + asample -= nsample; + if( app->getRngUniform() < asample ) ++nsample; + + EGS_Float aux1 = ct_max - ct_min, aux2 = 1 - beta*ct_max; + EGS_Float wt = pi.wt/nbrspl; + EGS_Float sinpsi, sindel, cosdel; + bool need_rotation; + sinpsi = u.x*u.x + u.y*u.y; + if( sinpsi > 1e-20 ) { + sinpsi = sqrt(sinpsi); + need_rotation = true; + cosdel = u.x/sinpsi; + sindel = u.y/sinpsi; + } else need_rotation = false; + int ip = np; + int ib = 0; + + if( w1 < w2 ) { + for(int j=0; jgetRngUniform(), cost; + if( app->getIbrdst() == 1 ) { + EGS_Float tmp = cmin*cmax/(rnno*cmin+(1-rnno)*cmax); + cost = (1 - tmp*tmp)/beta; + } + else { + rnno *= aux1; + cost = (ct_min*aux2+rnno)/(aux2+beta*rnno); + } + EGS_Float sint = 1 - cost*cost; + sint = sint > 0 ? sqrt(sint) : 0; + EGS_Float cphi,sphi; + app->getRngAzimuth(cphi,sphi); + EGS_Float un,vn,wn; + if( need_rotation ) { + EGS_Float us = sint*cphi, vs = sint*sphi; + un = u.z*cosdel*us - sindel*vs + u.x*cost; + vn = u.z*sindel*us + cosdel*vs + u.y*cost; + wn = u.z*cost - sinpsi*us; + } else { + un = sint*cphi; + vn = sint*sphi; + wn = u.z*cost; + } + int ns = 0; + if( wn > 0 ) { + EGS_Float aux = (ssd - x.z)/wn; + EGS_Float x1 = x.x + un*aux, y1 = x.y + vn*aux; + if( x1*x1 + y1*y1 < fs*fs ) ns = 1; + } + /* + if( !ns ) { + if( app->getRngUniform()*nbrspl < 1 ) + { + //odd: we allow the generation of a fat photon here + //this is a hangover from the beampp implementation...comment out for now + ns = nbrspl; + the_fat = true; + } + } + */ + if( ns > 0 ) { + if( ++ip >= app->getMxstack() ) egsFatal(dbs_err_msg, + "smartBrems",app->getMxstack()); + if( app->getIbrdst() == 1 ) + y2_KM[ib++] = beta*(1+beta)*(tau+1)*(tau+1)*(1-cost); + + EGS_Particle p; + p.x = x; + p.u = EGS_Vector(un,vn,wn); + p.ir = irl; + p.wt = wt*ns; + p.latch = latch; + if(the_fat) p.latch = p.latch | (1 << 0); + p.q = 0; + + particle_stack.emplace_back(p); + + /* need to figure this out + the_extra_stack->icreate[ip] = the_extra_stack->int_num; + the_extra_stack->pid[ip] = ++the_extra_stack->pidI; + the_extra_stack->iweight[ip] = ns; + */ + } + } + } + else { + for(int j=0; jgetPoint(fs,x1,y1,iw); + //else { + do { + x1 = 2*app->getRngUniform()-1; + y1 = 2*app->getRngUniform()-1; + } while ( x1*x1 + y1*y1 > 1 ); + x1 *= fs; + y1 *= fs; + iw = 1; + //} + EGS_Float un = x1 - x.x, vn = y1 - x.y, wn = d; + EGS_Float dist = sqrt(un*un + vn*vn + wn*wn); + EGS_Float disti = 1/dist; + un *= disti; + vn *= disti; + wn *= disti; + EGS_Float cost = u.x*un + u.y*vn + u.z*wn; + EGS_Float aux = dmin*disti; + EGS_Float rejf; + if( app->getIbrdst() == 1 ) { + rejf = aux2/(1-beta*cost); + rejf *= sqrt(rejf)*aux*aux*aux; + } + else { + rejf = aux2/(1-beta*cost)*aux; + rejf *= rejf*aux; + } + + if( app->getRngUniform() < rejf ) { + if( ++ip >= app->getMxstack() ) egsFatal(dbs_err_msg, + "smartBrems",app->getMxstack()); + if( app->getIbrdst() == 1 ) { + EGS_Float y2 = beta*(1+beta)*(tau+1)*(tau+1)*(1-cost); + if( y2 < 0 ) y2 = 0; + y2_KM[ib++] = y2; + } + + EGS_Particle p; + p.x = x; + p.u = EGS_Vector(un,vn,wn); + p.ir = irl; + p.wt = wt*iw; + p.latch = latch; + p.q = 0; + + particle_stack.emplace_back(p); + + //have to figure out what to do below + /* + the_extra_stack->icreate[ip] = the_extra_stack->int_num; + the_extra_stack->pid[ip] = ++the_extra_stack->pidI; + the_extra_stack->iweight[ip] = iw; + */ + } + } + } + + //below generates a single phat photon directed away from the splitting field + int ireal = -1; + bool take_it = true; + EGS_Float un=0,vn=0,wn=1; + if( app->getIbrdst() != 1 ) { + aux1 = 2; + aux2 = 1 - beta; + EGS_Float rnno = app->getRngUniform()*aux1; + EGS_Float cost = (rnno-aux2)/(aux2+beta*rnno); + EGS_Float sint = 1 - cost*cost; + sint = sint > 0 ? sqrt(sint) : 0; + EGS_Float cphi,sphi; + app->getRngAzimuth(cphi,sphi); + if( need_rotation ) { + EGS_Float us = sint*cphi, vs = sint*sphi; + un = u.z*cosdel*us - sindel*vs + u.x*cost; + vn = u.z*sindel*us + cosdel*vs + u.y*cost; + wn = u.z*cost - sinpsi*us; + } else { + un = sint*cphi; + vn = sint*sphi; + wn = u.z*cost; + } + if( wn > 0) { + EGS_Float t = (ssd-x.z)/wn; + EGS_Float x1 = x.x + un*t, y1 = x.y + vn*t; + if( x1*x1 + y1*y1 <= fs*fs ) take_it = false; + } + } + if( take_it ) { + //put a high-weight particle directed away from the splitting field + if( ++ip >= app->getMxstack() ) egsFatal(dbs_err_msg, + "smartBrems",app->getMxstack()); + ireal = ip-np-1; //index of this high weight particle in the local particle_stack + + EGS_Particle p; + p.x = x; + p.u = EGS_Vector(un,vn,wn); + p.ir = irl; + p.wt = wt*nbrspl; + //label the particle as phat + p.latch = latch | (1 << 0); + p.q = 0; + + particle_stack.emplace_back(p); + + } + + //get the energies of the brem photons...and the initiating e- + getBremsEnergies(); + + if( ip > np && app->getIbrdst() == 1 ) { + ib = 0; + EGS_Float E = tau+1; + EGS_Float unorm = log(0.5*(E+1)); + for(int j=0; jgetRM(); + EGS_Float Ep = E - k; + EGS_Float r = Ep/E; + EGS_Float delta = k/(2*E*Ep); + delta *= delta; + EGS_Float arg1 = 1 + r*r, arg2 = arg1 + 2*r; + if( j != ireal ) { + EGS_Float y2 = y2_KM[ib++]; + EGS_Float x = k/tau; + EGS_Float f1,f2; + if( x <= 0.5 ) { + f1 = f_KM_a[imed][j_KM].interpolate(x); + f2 = f_KM_a[imed][j_KM+1].interpolate(x); + } else { + EGS_Float u = log(E-k)/unorm; + f1 = f_KM_b[imed][j_KM].interpolate(u); + f2 = f_KM_b[imed][j_KM+1].interpolate(u); + } + EGS_Float f = f1*p_KM + f2*q_KM; + EGS_Float a = y2 + 1, a2 = a*a; + EGS_Float rejf = 16*y2*r/a2 - arg2 - + (arg1 - 4*y2/a2*r)*log(delta+zbr_KM[imed]/a2); + rejf /= (f*sqrt(a)*f_max_KM); + if( app->getRngUniform() > rejf ) { + particle_stack[j].wt = 0; + particle_stack[j].E = 0; + } + } + else { + // this is the single high-weight photon + // now that we know the photon energy, we must sample + // a direction for this photon from 2BS + EGS_Float arg3 = -log(delta+zbr_KM[imed]); + EGS_Float rejmax = arg1*arg3-arg2; + EGS_Float y2max = 2*beta*(1+beta)*E*E; + EGS_Float z2maxi = sqrt(1+y2max); + EGS_Float rejf,y2,rnno; + do { + rnno = app->getRngUniform(); + y2 = app->getRngUniform(); + EGS_Float aux3 = z2maxi/(y2+(1-y2)*z2maxi); + y2 = aux3*aux3-1; + rnno *= rejmax*aux3; + EGS_Float aux3_4 = aux3*aux3; + aux3_4 *= aux3_4; + EGS_Float y2tst1 = r*y2/aux3_4; + EGS_Float aux4 = 16*y2tst1-arg2, + aux5 = arg1-4*y2tst1; + if( rnno < aux4 + aux5*arg3 ) break; + EGS_Float aux2=log(aux3_4/(delta*aux3_4+zbr_KM[imed])); + rejf = aux4+aux5*aux2; + } while (rejf < rnno); + EGS_Float cost = 1 - 2*y2/y2max; + bool take_it; + EGS_Float un,vn,wn; + if( w2<=w1 || ((costct_max) && w1 0 ? sqrt(sint) : 0; + EGS_Float cphi,sphi; + app->getRngAzimuth(cphi,sphi); + if( need_rotation ) { + EGS_Float us = sint*cphi, vs = sint*sphi; + un = u.z*cosdel*us - sindel*vs + u.x*cost; + vn = u.z*sindel*us + cosdel*vs + u.y*cost; + wn = u.z*cost - sinpsi*us; + } else { + un = sint*cphi; + vn = sint*sphi; + wn = u.z*cost; + } + if( w2 <= w1 && wn > 0) { + EGS_Float t = (ssd-x.z)/wn; + EGS_Float x1 = x.x + un*t, y1 = x.y + vn*t; + if( x1*x1 + y1*y1 <= fs*fs ) take_it = false; + } + } else take_it = false; + if( take_it ) { + particle_stack[j].u = EGS_Vector(un,vn,wn); + } + else { + particle_stack[j].wt = 0; + particle_stack[j].E = 0; + } + } + } + } + //now add particles to the stack + for (int i=0; iaddParticleToStack(particle_stack[i],dneari); + } + + return 0; +} + + +void EGS_RadiativeSplitting::getCostMinMax(const EGS_Vector &xx, const EGS_Vector &uu, + EGS_Float &ro, EGS_Float &ct_min, EGS_Float &ct_max) { + double d = ssd - xx.z, d2 = d*d; + double xo = xx.x, yo = xx.y; + double u = uu.x, v = uu.y, w = uu.z; + double ro2 = xo*xo + yo*yo; + ro = sqrt(ro2); + double r = fs; + double r2 = r*r; + double st2 = u*u + v*v; + double st = sqrt(st2); + + // handle odd cases st=0 and/or ro=0 + if( st2 < 1e-10 ) { + double dmin = sqrt(d2 + (r-ro)*(r-ro)), + dmax = sqrt(d2 + (r+ro)*(r+ro)); + if( w > 0 ) { + ct_max = ro <= r ? 1 : d*w/dmin; + ct_min = d*w/dmax; + } + else { + ct_max = d*w/dmax; + ct_min = ro <= r ? -1 : d*w/dmin; + } + return; + } + if( ro2 < 1e-8 ) { + double aux = 1/sqrt(d2 + r*r); + ct_max = (d*w + r*st)*aux; + ct_min = (d*w - r*st)*aux; + if( w ) { + double x1 = xo + u*d/w, y1 = yo + v*d/w; + if( x1*x1 + y1*y1 <= r2 ) { + if( w > 0 ) ct_max = 1; + else ct_min = -1; + } + } + return; + } + + EGS_Float dmin = ro <= fs ? d : sqrt(d2 + (fs-ro)*(fs-ro)); + EGS_Float dmax = sqrt(d2 + (fs+ro)*(fs+ro)); + EGS_Float aux = w*d - u*xo - v*yo; + ct_max = aux + fs*st; + if( ct_max > 0 ) ct_max /= dmin; + else ct_max /= dmax; + if( ct_max > 1 ) ct_max = 1; + ct_min = aux - fs*st; + if( ct_min > 0 ) ct_min /= dmax; + else ct_min /= dmin; + if( ct_min < -1 ) ct_min = -1; +} + +void EGS_RadiativeSplitting::getBremsEnergies() { + + //get energies of brems photons and intiating e- + //copied from beamnrc.mortran, which is, in turn, copied from subroutine BREMS in egsnrc.mortran + EGS_Float eie, ekin, ese, esg, brmin, waux, r1,ajj,br,rnno06,rnno07,delta,phi1,phi2,rejf,x1,y1,aux; + int nsample,ip,j,jj,l,l1, imed; + EGS_Float peie,pese,pesg; + + int np = app->getNp(); + EGS_Particle pi = app->getParticleFromStack(np); + EGS_Float dneari = app->getDnear(np); + + int irl = pi.ir; + imed = app->getMedium(irl); + + peie = pi.E; //i.e. we have not added particles to the stack yet + eie = peie; + if (eie < 50.0) { + l=0; //index -1 compared to Mortran + } + else + { + l=2; //index -1 compared to Mortran + } + l1 = l+1; + ekin = peie - app->getRM(); + brmin = app->getAp(imed)/ekin; + waux = log(ekin) - log(app->getAp(imed)); + + if (app->getIbrnist() == 1) + { + ajj = 1 + (waux + log(app->getAp(imed)) - app->getNbLemin(imed))*app->getNbDlei(imed); + jj = ajj; + ajj = ajj - jj; + if ( jj > app->getMxbres()) { + jj = app->getMxbres(); + ajj = -1; + } + } + for (int ibr = 0; ibr < particle_stack.size(); ibr++) + { + if (app->getIbrnist() == 1) + { + if (ekin > app->getNbEmin(imed)) + { + r1 = app->getRngUniform(); + if (r1 < ajj) { + j = jj+1; + } + else + { + j = jj; + } + j=j-1; //index -1 compared to Mortran + //maybe should use egs++ alias_sample function here + int mxbrxs = app->getMxbrxs(); + EGS_Float* f1 = app->getNbXdata(j,imed); + EGS_Float* f2 = app->getNbFdata(j,imed); + EGS_Float* f3 = app->getNbWdata(j,imed); + int* f4 = app->getNbIdata(j,imed); + br = app->callAliasSample1(mxbrxs,f1,f2,f3,f4); + } + else + { + br = app->getRngUniform(); + } + esg = app->getAp(imed)*exp(br*waux); + pesg = esg; + pese = peie - pesg; + ese = pese; + } + else + { + do + { + // use BH cross-sections + rnno06 = app->getRngUniform(); + rnno07 = app->getRngUniform(); + br = brmin*exp(rnno06*waux); + esg = ekin*br; + pesg = esg; + pese = peie - pesg; + ese = pese; + delta = esg/eie/ese*app->getDelcm(imed); + aux = ese/eie; + if (delta < 1) + { + phi1 = app->getDl1(l,imed)+delta* + (app->getDl2(l,imed)+app->getDl3(l,imed)); + phi2 = app->getDl1(l1,imed)+delta* + (app->getDl2(l1,imed)+app->getDl3(l1,imed)); + } + else + { + phi1 = app->getDl4(l,imed)+app->getDl5(l,imed)* + log(delta + app->getDl6(l,imed)); + phi2 = phi1; + } + rejf = (1+aux*aux)*phi1 - 2*aux*phi2/3; + + } while (rnno07 >= rejf); + } + particle_stack[ibr].E = pesg; + } + //now (re)set energy of initiating e- + pi.E = pese; + app->updateParticleOnStack(np,pi,dneari); +} + +void EGS_RadiativeSplitting::killThePhotons(EGS_Float fs, EGS_Float ssd, int n_split, int npstart, int kill_electrons) { + //an adaptation of the Mortran subroutine kill_the_photons found in beamnrc.mortran, beampp.mortran + if (npstart > app->getNp()) return; + int i_playrr = 0; + int idbs = npstart; + int np; + EGS_Float dnear = app->getDnear(idbs); + do { + EGS_Particle p = app->getParticleFromStack(idbs); + inverseTransformP(p,T); + if (p.q == 0) + { + i_playrr = 0; + if (p.u.z < 0) + { + i_playrr = 1; + } + else + { + EGS_Float dist = (ssd - p.x.z)/p.u.z; + EGS_Float r2 = (p.x.x+dist*p.u.x)*(p.x.x+dist*p.u.x) + (p.x.y+dist*p.u.y)*(p.x.y+dist*p.u.y); + if (r2 > fs*fs) i_playrr = 1; + } + if (i_playrr) + { + if (app->getRngUniform()*n_split > 1) + { + //kill the particle + app->deleteParticleFromStack(idbs); + } + else + { + //keep the particle and increase weight + p.wt = p.wt*n_split; + //set bit 0 of latch to mark as phat + p.latch = p.latch | (1 << 0); + transformP(p,T); + app->updateParticleOnStack(idbs,p,dnear); + idbs++; + } + } + else + { + //keep the particle + idbs++; + } + } + else + { + //this is a charged particle + if (kill_electrons == 0) + { + //keep it + idbs++; + } + else + { + if (app->getRngUniform()*n_split > 1) + { + //kill it + app->deleteParticleFromStack(idbs); + } + else + { + //keep the particle, increase weight and label as phat + p.wt = p.wt*n_split; + p.latch = p.latch | (1 << 0); + transformP(p,T); + app->updateParticleOnStack(idbs,p,dnear); + idbs++; + } + } + } + np = app->getNp(); + } while (idbs <= np); + + return; +} + +void EGS_RadiativeSplitting::selectAzimuthalAngle(EGS_Float &cphi, EGS_Float &sphi) +{ + //function to randomly select an azimuthal angle evenly distributed on 4pi + //returns cos and sin of angle + //C++ version of $SELECT-AZIMUTHAL-ANGLE macro in egsnrc.macros + EGS_Float xphi,yphi,xphi2,yphi2,rhophi2; + do { + xphi = app->getRngUniform(); + xphi = 2*xphi - 1; + xphi2 = xphi*xphi; + yphi = app->getRngUniform(); + yphi2 = yphi*yphi; + rhophi2 = xphi2 + yphi2; + } while (rhophi2 > 1); + rhophi2 = 1/rhophi2; + cphi = (xphi2 - yphi2)*rhophi2; + sphi = 2*xphi*yphi*rhophi2; + return; +} + +void EGS_RadiativeSplitting::uniformPhotons(int nsample, int n_split, EGS_Float fs, EGS_Float ssd, EGS_Float energy) +{ + //function to generate photons radiating isotropically from a point that are directed into the splitting + //field defined by fs, ssd + //modeled after the Mortran routine used in beamnrc, beampp + + //get properties of interacting particle (annihilating positron or radiative photon being split) + EGS_Particle pi = app->getParticleFromStack(app->getNp()); + inverseTransformP(pi,T); + EGS_Float x = pi.x.x; + EGS_Float y = pi.x.y; + EGS_Float z = pi.x.z; + EGS_Float dnear = app->getDnear(app->Np); + EGS_Float weight = pi.wt/n_split; + + //calculate min/max polar angles subtended by the splitting field + EGS_Float ro = sqrt(x*x+y*y); + EGS_Float d = ssd - z; + EGS_Float aux = (ro + fs)/d; + EGS_Float ct_min = 1/sqrt(1+aux*aux); + EGS_Float ct_max; + if (ro <= fs) + { + ct_max = 1; + } + else + { + aux = (fs - ro)/d; + ct_max = 1/sqrt(1+aux*aux); + } + + //calculate fraction of nsample photons that will have polar angles between ct_min/ct_max + EGS_Float an_split = 0.5*(ct_max - ct_min)*nsample; + int num_split = an_split; + an_split = an_split - num_split; + if (app->getRngUniform() < an_split) num_split += 1; + + //now create the num_split photons + //note: not all photons will go towards the circle and we will have to play russian roulette with these + EGS_Float cost, sint, cphi, sphi; + int ns,ip; + //delete interacting particle on stack so that we overwrite it + app->deleteParticleFromStack(app->getNp()); + //and update np + int np = app->getNp(); + for (int i=0; i< num_split; i++) + { + EGS_Particle p; + cost = ct_min + app->getRngUniform()*(ct_max-ct_min); + sint = 1-cost*cost; + if (sint > 0) + { + sint = sqrt(sint); + } + else + { + sint = 0; + } + selectAzimuthalAngle(cphi,sphi); + aux = d/cost*sint; + EGS_Float x_tmp = x + aux*cphi; + EGS_Float y_tmp = y + aux*sphi; + ns = 0; + if (x_tmp*x_tmp + y_tmp*y_tmp < fs*fs) + { + ns = 1; + } + else + { + //this is directed outside the field, play RR in next if block + ns = nsplit; + } + if (ns > 1) + { + if(app->getRngUniform()*ns>1) + { + ns=0; + } + } + if (ns > 0) + { + np++; + p.latch = pi.latch; + if (ns > 1) + { + // going outside field, label as phat + p.latch = p.latch | (1 << 0); + } + p.ir = pi.ir; + p.x = pi.x; + //stuff that we do not inherit + p.q = 0; + p.E = energy; + p.wt = weight*ns; + p.u = EGS_Vector(sint*cphi, sint*sphi, cost); + + transformP(p,T); + + app->addParticleToStack(p,dnear); + } + } + + //now generate (1-(ct_max-ct_min)/2*nsample photons directed outside the field + //model this by generating nsample/num_split photons over all polar angles and + //rejecting those with angles between ct_min and ct_max + num_split = nsample/n_split; + for (int i= 0; igetRngUniform()-1; + if (cost < ct_min || cost > ct_max) + { + sint = 1 - cost*cost; + if (sint > 0 ) + { + sint = sqrt(sint); + } + else + { + sint = 0; + } + selectAzimuthalAngle(cphi,sphi); + np++; + EGS_Particle p; + //label particle as phat + p.latch = pi.latch | (1 << 0); + p.ir = pi.ir; + p.x = pi.x; + //stuff that we do not inherit + p.q = 0; + p.E = energy; + p.wt = weight*n_split; + p.u = EGS_Vector(sint*cphi, sint*sphi, cost); + + transformP(p,T); + + app->addParticleToStack(p,dnear); + } + } + return; +} + +void EGS_RadiativeSplitting::doSmartCompton(int nint) +{ + //This is based on (i.e. copied from) the BEAMnrc implementation + + //get properties of interacting particle + int np = app->getNp(); + EGS_Particle pi = app->getParticleFromStack(np); + inverseTransformP(pi,T); + EGS_Float E = pi.E; + EGS_Vector x = pi.x; + EGS_Vector u = pi.u; + int irl=pi.ir, latch=pi.latch; + imed = app->getMedium(irl); + EGS_Float AP = app->getAp(imed); + EGS_Float dnear = app->getDnear(np); + //reduce weight of split particles + EGS_Float wt = pi.wt/nint; + + EGS_Float ct_min,ct_max,ro; + EGS_Float d = ssd - x.z; + getCostMinMax(x,u,ro,ct_min,ct_max); + + EGS_Float Ko = E/app->getRM(); + EGS_Float broi = 1+2*Ko, Ko2 = Ko*Ko; + EGS_Float alpha1_t = log(broi); + EGS_Float eps1_t = 1./broi, eps2_t = 1; + EGS_Float w2 = alpha1_t*(Ko2-2*Ko-2)+(eps2_t-eps1_t)* + (1./eps1_t/eps2_t + broi + Ko2*(eps1_t+eps2_t)/2); + EGS_Float eps12_t = eps1_t*eps1_t; + EGS_Float alpha2_t = (eps2_t*eps2_t-eps12_t); + EGS_Float alpha_t = alpha1_t/(alpha1_t+alpha2_t/2); + EGS_Float eps1 = 1./(1+Ko*(1-ct_min)), eps2 = 1./(1+Ko*(1-ct_max)); + EGS_Float eps1_0 = eps1, eps2_0 = eps2; + EGS_Float alpha1 = log(eps2/eps1); + EGS_Float w1 = alpha1*(Ko2-2*Ko-2)+(eps2-eps1)*(1./eps1/eps2 + broi + + Ko2*(eps1+eps2)/2); + EGS_Float eps12 = eps1*eps1, alpha2 = (eps2*eps2-eps12); + EGS_Float alpha = alpha1/(alpha1+alpha2/2); + EGS_Float rej1 = 1-(1-eps1)*(broi*eps1-1)/(Ko*Ko*eps1*(1+eps1*eps1)); + EGS_Float rej2 = 1-(1-eps2)*(broi*eps2-1)/(Ko*Ko*eps2*(1+eps2*eps2)); + EGS_Float rejmax = max(rej1,rej2); + + //determine no. of interactions to sample + EGS_Float wc = w1/w2; + EGS_Float asample = wc*nint; + int nsample = (int) asample; + asample -= nsample; + if( app->getRngUniform() < asample ) ++nsample; + + // prepare rotations--not totally sure why this is needed + EGS_Float sinpsi, sindel, cosdel; + bool need_rotation; + sinpsi = u.x*u.x + u.y*u.y; + if( sinpsi > 1e-20 ) { + sinpsi = sqrt(sinpsi); + need_rotation = true; + cosdel = u.x/sinpsi; + sindel = u.y/sinpsi; + } else need_rotation = false; + + // + // sample interactions towards circle + // + + EGS_Float br,sint,cost,cphi,sphi; //declared out here because they are also used for the electron below + EGS_Particle p; + + for(int j=0; jgetRngUniform() < alpha ) + br = eps1*exp(alpha1*app->getRngUniform()); + else + br = sqrt(eps12 + app->getRngUniform()*alpha2); + temp = (1-br)/(Ko*br); + sint = temp*(2-temp); + rejf = 1 - br*sint/(1+br*br); + } while ( app->getRngUniform()*rejmax > rejf ); + + if ( temp < 2 ) + { + cost = 1 - temp; + sint = sqrt(sint); + } + else + { + cost = -1; + sint = 0; + } + app->getRngAzimuth(cphi,sphi); + + int ns=0; + if (j==nsample-1) + { + //potential phat photon directed away from the field + if (br > eps1_0 && br < eps2_0) + { + break; // last time through the loop anyway, but do not generate this phat photon + } + ns = nint; + } + + EGS_Float un,vn,wn; + if( need_rotation ) + { + EGS_Float us = sint*cphi, vs = sint*sphi; + un = u.z*cosdel*us - sindel*vs + u.x*cost; + vn = u.z*sindel*us + cosdel*vs + u.y*cost; + wn = u.z*cost - sinpsi*us; + } + else + { + un = sint*cphi; + vn = sint*sphi; + wn = u.z*cost; + } + + if (j 0 ) { + EGS_Float aux = (ssd - x.z)/wn; + EGS_Float x1 = x.x + un*aux, y1 = x.y + vn*aux; + if( x1*x1 + y1*y1 < fs*fs ) ns = 1; + } + if (ns > 1) + { + //potentially a 2nd phat photon? + if (app->getRngUniform()*ns > 1) + { + ns = 0; + } + } + } + + if( ns > 0 ) { + //add the photon to the stack + p.x = x; + p.u = EGS_Vector(un,vn,wn); + p.ir = irl; + p.wt = wt*ns; + p.latch = latch; + if (ns > 1) + { + //label photon as phat + p.latch = latch | (1 << 0); + } + p.q = 0; + p.E = E*br; + + transformP(p,T); + + app->addParticleToStack(p,dnear); + + } + } + + //now the electron--Note: br, cost will have the values set for the phat photon + EGS_Float Eelec = E*(1-br); + EGS_Float aux = 1 + br*br - 2*br*cost; + EGS_Float un=u.x,vn=u.y,wn=u.z; + if( aux > 1e-8 ) { + cost = (1-br*cost)/sqrt(aux); + sint = 1-cost*cost; + if( sint > 0 ) sint = -sqrt(sint); + else sint = 0; + EGS_Float us = sint*cphi, vs = sint*sphi; + un = u.z*cosdel*us - sindel*vs + u.x*cost; + vn = u.z*sindel*us + cosdel*vs + u.y*cost; + wn = u.z*cost - sinpsi*us; + } + p.x = x; + p.u = EGS_Vector(un,vn,wn); + p.ir = irl; + p.wt = wt*nint; + p.latch = latch | (1 << 0); + p.q = -1; + p.E = Eelec + app->getRM(); + //replace original interacting photon with the e- + + transformP(p,T); + + app->updateParticleOnStack(np,p,dnear); +} + +void EGS_RadiativeSplitting::initSmartKM(EGS_Float Emax) { + + int nmed = app->getNmed(); + if( nmed < 1 ) return; + EGS_Float logEmax = log(Emax); + + y2_KM = new EGS_Float [nsplit]; + f_KM_a = new EGS_Interpolator* [nmed]; + f_KM_b = new EGS_Interpolator* [nmed]; + f_KM_max = new EGS_Interpolator [nmed]; + zbr_KM = new EGS_Float [nmed]; + a_KM = new EGS_Float [nmed]; + b_KM = new EGS_Float [nmed]; + + DBS_Aux the_aux; + the_aux.np = 128; + EGS_Float b = 0.5; + int nx = 257; + EGS_Float ddx = 0.00390625; + for(int imed=0; imedgetAp(imed); + EGS_Float logAP = log(AP); + int ne = 1 + (int) ((logEmax-logAP)/0.12); + a_KM[imed] = ((EGS_Float)ne-1)/(logEmax-logAP); + b_KM[imed] = -a_KM[imed]*logAP; + f_KM_a[imed] = new EGS_Interpolator [ne]; + f_KM_b[imed] = new EGS_Interpolator [ne]; + EGS_Float *tmp_fmax = new EGS_Float [ne]; + EGS_Float dle = (logEmax-logAP)/(ne-1); + int imedm = imed+1; + zbr_KM[imed] = app->getZbrang(imed); + the_aux.Z13 = zbr_KM[imed]; + for(int ie=0; iegetRM() + 1; + the_aux.E = E; + the_aux.le = log(E - b*(E-1)); + f_KM_a[imed][ie].initialize(0,b,dbs_func_KM,&the_aux,1025,9e-4); + f_KM_b[imed][ie].initialize(0,1,dbs_func_KM1,&the_aux,1025,9e-4); + EGS_Float fmax = 0; + for(int j=0; j fmax ) fmax = f; + } + tmp_fmax[ie] = fmax; + } + f_KM_max[imed].initialize(ne,logAP,logEmax,tmp_fmax); + //egsWarning("---deleting in initSmartKM()\n"); + delete [] tmp_fmax; + } + nmed_KM = nmed; + if( nsplit > 0 ) y2_KM = new EGS_Float [nsplit]; +} + +void EGS_RadiativeSplitting::doUphi21(EGS_Float sinthe, EGS_Float costhe, EGS_Vector &u) +{ + //This is equivalent to calling the Mortran routine UPHI(2,1) and is used + //to adjust particle angles after a Rayleigh event + EGS_Float cosphi,sinphi; + selectAzimuthalAngle(cosphi,sinphi); + EGS_Float a=u.x; + EGS_Float b=u.y; + EGS_Float c=u.z; + EGS_Float sinps2 = a*a+b*b; + if (sinps2 < 1e-20) + { + //small polar angle case + u.x = sinthe*cosphi; + u.y = sinthe*sinphi; + u.z = c*costhe; + } + else + { + EGS_Float sinpsi = sqrt(sinps2); + EGS_Float us = sinthe*cosphi; + EGS_Float vs = sinthe*sinphi; + EGS_Float sindel = b/sinpsi; + EGS_Float cosdel = a/sinpsi; + u.x = c*cosdel*us - sindel*vs + a*costhe; + u.y = c*sindel*us + cosdel*vs + b*costhe; + u.z = -sinpsi*us + c*costhe; + } + return; +} + //********************************************************************* // Process input for this ausgab object // @@ -103,14 +1672,89 @@ extern "C" { return 0; } + //get type of splitting + vector allowed_split_type; + allowed_split_type.push_back("uniform"); + allowed_split_type.push_back("directional"); + //TODO: Implement beampp-style DBS + int split_type; + string str; + if(input->getInput("splitting type",str) < 0) + { + split_type = EGS_RadiativeSplitting::URS; + } + else + { + split_type = input->getInput("splitting type",allowed_split_type,-1); + if ( split_type < EGS_RadiativeSplitting::URS || + split_type > EGS_RadiativeSplitting::DRSf ) + { + egsFatal("\nEGS_RadiativeSplitting: Invalid splitting type.\n"); + } + } EGS_Float nsplit = 1.0; + EGS_Float fs,ssd,zrr=0; + int irad = 0; + vector esplit_reg; + EGS_AffineTransform *t; int err = input->getInput("splitting",nsplit); + if (err) + { + egsWarning("\nEGS_RadiativeSplitting: Invalid splitting no. Radiative splitting will be turned off.\n"); + } + if ( split_type == EGS_RadiativeSplitting::DRS || + split_type == EGS_RadiativeSplitting::DRSf ) + { + int err01 = input->getInput("field size",fs); + if (err01) + { + egsFatal("\nEGS_RadiativeSplitting: Missing/invalid input for splitting field radius.\n"); + } + int err02 = input->getInput("ssd",ssd); + if (err02) + { + egsFatal("\nEGS_RadiativeSplitting: Missing/invalid input for splitting field ssd.\n"); + } + int err03 = input->getInput("e-/e+ split regions",esplit_reg); + if (err03) + { + egsWarning("\nEGS_RadiativeSplitting: Missing/invalid input for e+/e- splitting region.\nCharged particles will not be split.\n"); + } + else + { + if(!input->getInput("radially redistribute e-/e+",str)) + { + vector allowed_irad; + allowed_irad.push_back("no"); + allowed_irad.push_back("yes"); + irad = input->getInput("radially redistribute e-/e+",allowed_irad,-1); + if (irad < 0) + { + egsWarning("\nEGS_RadiativeSplitting: Invalid input for e-/e+ radial redistribution. Will not radially redistribute split charged particles.\n"); + irad = 0; + } + } + int err04 = input->getInput("Z of russian roulette plane",zrr); + if (err04) + { + egsWarning("\nEGS_RadiativeSplitting: Missing/invalid input for Z of Russian Roulette plane. Defaults to 0.\n"); + } + } - //================================================= + //option to transform DBS splitting cone + t = EGS_AffineTransform::getTransformation(input); + } + //================================================= /* Setup radiative splitting object with input parameters */ EGS_RadiativeSplitting *result = new EGS_RadiativeSplitting("",f); + result->setSplitType(split_type); result->setSplitting(nsplit); + if (split_type == EGS_RadiativeSplitting::DRS || + split_type == EGS_RadiativeSplitting::DRSf) { + result->initDBS(fs,ssd,esplit_reg,irad,zrr,t); + delete t; + } result->setName(input); return result; } 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..0a3c8f0d2 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 @@ -23,7 +23,7 @@ # # Author: Ernesto Mainegra-Hing, 2018 # -# Contributors: +# Contributors: Blake Walters, 2021 # ############################################################################### # @@ -31,7 +31,9 @@ # # TODO: # -# - Add directional radiative splitting (DRS) +# - Testing/debugging directional radiative splitting (DRS) +# - Implement e- splitting +# - Implement beampp-style DBS (DRSf) # ############################################################################### */ @@ -49,25 +51,28 @@ #include "egs_application.h" #include "egs_scoring.h" #include "egs_base_geometry.h" +#include "egs_rndm.h" +#include "egs_interpolator.h" +#include "egs_transformations.h" #ifdef WIN32 - #ifdef BUILD_RADIATIVE_SPLITTING_DLL - #define EGS_RADIATIVE_SPLITTING_EXPORT __declspec(dllexport) - #else - #define EGS_RADIATIVE_SPLITTING_EXPORT __declspec(dllimport) - #endif - #define EGS_RADIATIVE_SPLITTING_LOCAL +#ifdef BUILD_RADIATIVE_SPLITTING_DLL +#define EGS_RADIATIVE_SPLITTING_EXPORT __declspec(dllexport) +#else +#define EGS_RADIATIVE_SPLITTING_EXPORT __declspec(dllimport) +#endif +#define EGS_RADIATIVE_SPLITTING_LOCAL #else - #ifdef HAVE_VISIBILITY - #define EGS_RADIATIVE_SPLITTING_EXPORT __attribute__ ((visibility ("default"))) - #define EGS_RADIATIVE_SPLITTING_LOCAL __attribute__ ((visibility ("hidden"))) - #else - #define EGS_RADIATIVE_SPLITTING_EXPORT - #define EGS_RADIATIVE_SPLITTING_LOCAL - #endif +#ifdef HAVE_VISIBILITY +#define EGS_RADIATIVE_SPLITTING_EXPORT __attribute__ ((visibility ("default"))) +#define EGS_RADIATIVE_SPLITTING_LOCAL __attribute__ ((visibility ("hidden"))) +#else +#define EGS_RADIATIVE_SPLITTING_EXPORT +#define EGS_RADIATIVE_SPLITTING_LOCAL +#endif #endif @@ -81,12 +86,29 @@ events. This ausgab object is specified via: :start ausgab object: library = egs_radiative_splitting name = some_name - splitting = n_split + splitting type = uniform (default) or directional (based on BEAMnrc DBS) + splitting = the splitting number (n_split) + + The following inputs apply to directional splitting only: + field size = radius of splitting field centred on Z-axis (cm) -- required + ssd = distance of splitting field from Z=0 along +Z-axis (cm) -- required + e-/e+ split regions = region number(s) for e-/e+ splitting. On entering this(ese) region(s), phat charged particles will be split n_split times. + If set to 0 or omitted, charged particles will not be split. + radially redistribute e-/e+ = "yes" or "no" (default) -- if "yes", evenly distribute split e-/e+ in a circle of radius sqrt(x(np)^2+y(np)^2) about the Z-axis + Z of russian roulette plane = Z below which russian roulette is not played on low-weight charged particles resulting from e-/e+ splitting (cm) -- optimally set to some value < Z at which particles enter the splitting region(s) + Below is an optional input for an affine transform to be applied to the directional splitting cone defined by field size and ssd: + :start transformation: + rotation = optional 3D rotation vector to be applied to the splitting cone (applied first) + translation = optional 3D translation to be applied to the splitting cone + :stop transformation: + :stop ausgab object: \endverbatim TODO: - - Add directional radiative splitting (DRS) + - Testing/debugging directional radiative splitting (DRS) + - Implement e-/e+ splitting + - Implement beampp-style DRS (DRSf) */ @@ -95,10 +117,10 @@ class EGS_RADIATIVE_SPLITTING_EXPORT EGS_RadiativeSplitting : public EGS_AusgabO public: /*! Splitting algortihm type */ - enum Type { + enum { URS, // EGSnrc Uniform Radiative Splitting - DRS, // Directional Radiative Splitting - DRSf // Directional Radiative Splitting (BEAMnrc) + DRS, // Directional Radiative Splitting (based on BEAMnrc DBS) + DRSf // Directional Radiative Splitting (based on beampp DBS) }; EGS_RadiativeSplitting(const string &Name="", EGS_ObjectFactory *f = 0); @@ -107,21 +129,147 @@ class EGS_RADIATIVE_SPLITTING_EXPORT EGS_RadiativeSplitting : public EGS_AusgabO void setApplication(EGS_Application *App); + void initializeData() { + //set bit 0 of the first particle in each history to 1 + //Until a more general solution is implemented (extra_stack), this + //bit determines whether the particle (and its descendents) has(ve) + //been split (i.e. is/are phat) or not + //Note: This means we lose the ability of this bit as a flag for + //brem events + EGS_Particle p = app->getParticleFromStack(0); + int latch = p.latch; + latch = latch | (1 << 0); + app->setLatch(latch); + } + + int doInteractions(int iarg, int &killed); + + int doSmartBrems(); + + void doSmartCompton(int nsample); + + void getCostMinMax(const EGS_Vector &xx, const EGS_Vector &uu, + EGS_Float &ro, EGS_Float &ct_min, EGS_Float &ct_max); + + void getBremsEnergies(); + + void killThePhotons(EGS_Float fs, EGS_Float ssd, int n_split, int npstart, int kill_electrons); + + void selectAzimuthalAngle(EGS_Float &cphi, EGS_Float &sphi); + + void uniformPhotons(int nsample, int n_split, EGS_Float fs, EGS_Float ssd, EGS_Float energy); + + void initSmartKM(EGS_Float Emax); + + void doUphi21(EGS_Float sinthe, EGS_Float costhe, EGS_Vector &u); + void setSplitting(const int &n_s) { nsplit = n_s; }; + void setSplitType(const int &type) { + split_type = type; + }; + + void initDBS(const EGS_Float field_rad, const EGS_Float field_ssd, const vector splitreg, const int irad, const EGS_Float zrr, const EGS_AffineTransform *t); + + bool needsCall(EGS_Application::AusgabCall iarg) const override { + if ( split_type == DRS || split_type == DRSf ) { + if ( iarg == EGS_Application::BeforeBrems || + iarg == EGS_Application::BeforeAnnihFlight || + iarg == EGS_Application::BeforeAnnihRest || + iarg == EGS_Application::BeforePair || + iarg == EGS_Application::BeforeCompton || + iarg == EGS_Application::BeforePhoto || + iarg == EGS_Application::BeforeRayleigh || + iarg == EGS_Application::FluorescentEvent || + iarg == EGS_Application::BeforeTransport) { + return true; + } + if ( ireg_esplit.size()>0 && iarg == EGS_Application::AfterTransport) { + return true; + } + } + return false; + }; + int processEvent(EGS_Application::AusgabCall iarg) { + if ( split_type > URS && iarg >= EGS_Application::BeforeTransport ) + { + if( !doInteractions(iarg,killed) ) + { + return 0; + } + } return 0; }; + int processEvent(EGS_Application::AusgabCall iarg, int ir) { + if (split_type > URS && iarg >= EGS_Application::BeforeTransport) + { + if( !doInteractions(iarg,killed) ) + { + return 0; + } + } return 0; }; + //function to rotate/translate particle position and rotate particle direction + void transformP(EGS_Particle &p, const EGS_AffineTransform *t) { + if (t) + { + t->transform(p.x); + t->rotate(p.u); + } + } + + //function to inverse rotate/translate particle position and inverse rotate direction + void inverseTransformP(EGS_Particle &p, const EGS_AffineTransform *t) { + if (t) + { + t->inverseTransform(p.x); + t->rotateInverse(p.u); + } + } + protected: + + int split_type; //0 = uniform, 1 = DBS, 2 = BEAMnrc DBS /* Maximum splitting limited to 2,147,483,647 */ int nsplit; + int killed; //no of photons killed + EGS_Float be_factor = 0; //Brems enhancement factor--currently set to zero and not used + EGS_Float fs; //radius of splitting field + EGS_Float ssd; //ssd at which splitting field is defined + vector ireg_esplit; //numbers of regions on entering which phat charged particles are split + int irad_esplit; //set to 1 to radially redistribute split e-/e+ + EGS_Float zrr_esplit; //Z value below which Russian Roulette will not be played with split e-/e+ + + bool use_cyl_sym = false; //set to true to use cylindrical symmetry, hard coded as false for now + EGS_Float zcyls; //Z below which cylindrical symmetry does not exist + + //interpolators for smart brems estimates of KM distribution + EGS_Interpolator *f_KM_max; + EGS_Interpolator **f_KM_a; + EGS_Interpolator **f_KM_b; + EGS_Float *a_KM, *b_KM; + EGS_Float *y2_KM; + EGS_Float *zbr_KM; + int nmed_KM; + + bool flag_fluor = false; + + EGS_AffineTransform *T; //optional transformation of the splitting cone + + vector particle_stack; //store a stack of brems particles in do_smart_brems + vector dnear_stack; //similar for dnear + + int imed; + const char *dbs_err_msg = + "Stack size exceeded in BEAMpp_DBS::%s()\n" + "Increase MXSTACK (currently %d) in array_sizes.h and retry\n"; }; #endif diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/egs_radiative_splitting.macros b/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/egs_radiative_splitting.macros new file mode 100644 index 000000000..47d0d82c8 --- /dev/null +++ b/HEN_HOUSE/egs++/ausgab_objects/egs_radiative_splitting/egs_radiative_splitting.macros @@ -0,0 +1,42 @@ +%C80 +"#############################################################################" +" " +" EGSnrc common block reorganization for C interface " +" Copyright (C) 2015 National Research Council Canada " +" " +" This file is part of EGSnrc. " +" " +" EGSnrc 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. " +" " +" EGSnrc 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. " +" " +" You should have received a copy of the GNU Affero General Public License " +" along with EGSnrc. If not, see . " +" " +"#############################################################################" +" " +" Author: Blake Walters 2021 " +" " +"#############################################################################" +" " +" This is a set of macros used to redefine some macros in " +" egs_c_interface2.macros as ; (i.e. null). This is coding that is not " +" required for the subset of egsnrc.mortran used by egs_radiative_splitting " +" and setting these blocks to null means that we do not have redefine C++ " +" versions of some subroutines " +"#############################################################################" + +REPLACE {$CALL-HOWNEAR(#);} WITH {;} + +REPLACE {$CALL-HOWFAR-IN-ELECTR;} WITH {;} + +REPLACE {$CALL-HOWFAR-IN-PHOTON;} WITH {;} + +REPLACE {$start_new_particle;} WITH {;} + diff --git a/HEN_HOUSE/egs++/egs_advanced_application.cpp b/HEN_HOUSE/egs++/egs_advanced_application.cpp index 5b55abf4f..3443c9254 100644 --- a/HEN_HOUSE/egs++/egs_advanced_application.cpp +++ b/HEN_HOUSE/egs++/egs_advanced_application.cpp @@ -210,6 +210,16 @@ int EGS_AdvancedApplication::initEGSnrcBackEnd() { return 0; } +//EGSnrc functions needed for radiative splitting ausgab object +extern __extc__ void F77_OBJ_(brems,BREMS)(); +extern __extc__ void F77_OBJ_(annih,ANNIH)(); +extern __extc__ void F77_OBJ_(annih_at_rest,ANNIH_AT_REST)(); +extern __extc__ void F77_OBJ_(photo,PHOTO)(); +extern __extc__ void F77_OBJ(pair,PAIR)(); +extern __extc__ void F77_OBJ_(compt,COMPT)(); +extern __extc__ void F77_OBJ_(egs_rayleigh_sampling,EGS_RAYLEIGH_SAMPLING)(int *medium, EGS_Float *e, EGS_Float *gle, EGS_I32 *lgle, EGS_Float *costhe, EGS_Float *sinthe); +extern __extc__ EGS_Float F77_OBJ_(alias_sample1,ALIAS_SAMPLE1)(int &mxbrxs, EGS_Float* nb_xdata, EGS_Float* nb_fdata, EGS_Float* nb_wdata, int* nb_idata); + /* The following 2 functions were used in an attempt to redirect all fortran I/O to use egsInformationn a degsWarning. Unfortunately, there were segmentation violations in the @@ -859,11 +869,11 @@ void EGS_AdvancedApplication::setEIIData(EGS_I32 len) { } #ifdef GDEBUG - #define MAX_STEP 100000 - EGS_Vector steps_x[MAX_STEP], steps_u[MAX_STEP]; - int steps_ireg[MAX_STEP], steps_inew[MAX_STEP]; - EGS_Float steps_ustepi[MAX_STEP], steps_ustepf[MAX_STEP]; - int steps_n = 0; +#define MAX_STEP 100000 +EGS_Vector steps_x[MAX_STEP], steps_u[MAX_STEP]; +int steps_ireg[MAX_STEP], steps_inew[MAX_STEP]; +EGS_Float steps_ustepi[MAX_STEP], steps_ustepf[MAX_STEP]; +int steps_n = 0; #endif int EGS_AdvancedApplication::shower() { @@ -886,6 +896,9 @@ int EGS_AdvancedApplication::shower() { #ifdef GDEBUG steps_n = 0; #endif + for (int j=0; jinitializeData(); + } egsShower(); return 0; } @@ -1261,6 +1274,260 @@ void EGS_AdvancedApplication::setLatch(int latch) { the_stack->latch[np] = latch; } +//************************************************************ +// Utility functions for ausgab egs_radiative_splitting objects +//************************************************************ + +//necessary RNG functions +EGS_Float EGS_AdvancedApplication::getRngUniform() { + return rndm->getUniform(); +} + +void EGS_AdvancedApplication::getRngAzimuth(EGS_Float &cphi, EGS_Float &sphi) { + rndm->getAzimuth(cphi,sphi); +} + +//add a particle (+ dnear) to the top of the stack and increment np +void EGS_AdvancedApplication::addParticleToStack(EGS_Particle p, EGS_Float dnear) { + int np = the_stack->np; //np+1 + the_stack->E[np] = p.E; + the_stack->wt[np] = p.wt; + the_stack->x[np] = p.x.x; + the_stack->y[np] = p.x.y; + the_stack->z[np] = p.x.z; + the_stack->u[np] = p.u.x; + the_stack->v[np] = p.u.y; + the_stack->w[np] = p.u.z; + the_stack->iq[np] = p.q; + the_stack->ir[np] = p.ir+2; + the_stack->latch[np] = p.latch; + the_stack->dnear[np] = dnear; + the_stack->np++; +} + +//get dnear from position np in the stack +EGS_Float EGS_AdvancedApplication::getDnear(int np) { + return the_stack->dnear[np]; +} + +//set npold of stack +void EGS_AdvancedApplication::setNpold(int npold) { + the_stack->npold = npold+1; +} + +//get total no. of media in simulation +int EGS_AdvancedApplication::getNmed() { + return nmed; +} + +//delete particle at stack position ip +//replace with data at position np and reduce np by 1 +void EGS_AdvancedApplication::deleteParticleFromStack(int ip) { + int np = the_stack->np-1; + if (ip < np) + { + the_stack->iq[ip] = the_stack->iq[np]; + the_stack->E[ip] = the_stack->E[np]; + the_stack->x[ip] = the_stack->x[np]; + the_stack->y[ip] = the_stack->y[np]; + the_stack->z[ip] = the_stack->z[np]; + the_stack->u[ip] = the_stack->u[np]; + the_stack->v[ip] = the_stack->v[np]; + the_stack->w[ip] = the_stack->w[np]; + the_stack->wt[ip] = the_stack->wt[np]; + the_stack->ir[ip] = the_stack->ir[np]; + the_stack->latch[ip] = the_stack->latch[np]; + the_stack->dnear[ip] = the_stack->dnear[np]; + } + the_stack->np--; + return; +} + +//retrieve particle information at stack position ip +EGS_Particle EGS_AdvancedApplication::getParticleFromStack(int ip) { + EGS_Particle p; + p.q = the_stack->iq[ip]; + p.E = the_stack->E[ip]; + p.latch = the_stack->latch[ip]; + p.ir = the_stack->ir[ip]-2; + p.wt = the_stack->wt[ip]; + p.x = EGS_Vector(the_stack->x[ip],the_stack->y[ip],the_stack->z[ip]); + p.u = EGS_Vector(the_stack->u[ip],the_stack->v[ip],the_stack->w[ip]); + return p; +} + +//update particle info (+ dnear) at stack position ip +void EGS_AdvancedApplication::updateParticleOnStack(int ip, EGS_Particle p, EGS_Float dnear) { + the_stack->iq[ip] = p.q; + the_stack->E[ip] = p.E; + the_stack->latch[ip] = p.latch; + the_stack->ir[ip] = p.ir+2; + the_stack->wt[ip] = p.wt; + the_stack->x[ip] = p.x.x; + the_stack->y[ip] = p.x.y; + the_stack->z[ip] = p.x.z; + the_stack->u[ip] = p.u.x; + the_stack->v[ip] = p.u.y; + the_stack->w[ip] = p.u.z; + the_stack->dnear[ip] = dnear; + return; +} + +//return the value of gle +EGS_Float EGS_AdvancedApplication::getGle() { + return the_epcont->gle; +} + +//return lgle +int EGS_AdvancedApplication::getLgle(EGS_Float gle, int med) { + return the_photin->ge1[med]*gle+the_photin->ge0[med]; +} + +//return the value of the_xoptions->ibrdst +int EGS_AdvancedApplication::getIbrdst() { + return the_xoptions->ibrdst; +} + +//return the value of the_xoptions->ibcmp +int EGS_AdvancedApplication::getIbcmp() { + return the_xoptions->ibcmp; +} + +//return the value of the_thresh->ap[imed] +EGS_Float EGS_AdvancedApplication::getAp(int imed) { + return the_thresh->ap[imed]; +} + +//return the value of the_xoptions->ibr_nist +int EGS_AdvancedApplication::getIbrnist() { + return the_xoptions->ibr_nist; +} + +//return the value of the_nist_brems->nb_lemin[imed] +EGS_Float EGS_AdvancedApplication::getNbLemin(int imed) { + return the_nist_brems->nb_lemin[imed]; +} + +//return the value of the_nist_brems->nb_dlei[imed] +EGS_Float EGS_AdvancedApplication::getNbDlei(int imed) { + return the_nist_brems->nb_dlei[imed]; +} + +//return the value of the_nist_brems->nb_emin[imed] +EGS_Float EGS_AdvancedApplication::getNbEmin(int imed) { + return the_nist_brems->nb_emin[imed]; +} + +//return the_nist_brems->nb_xdata[j,imed] +EGS_Float* EGS_AdvancedApplication::getNbXdata(int j, int imed) { + return the_nist_brems->nb_xdata[imed][j]; +} + +//return the_nist_brems->nb_fdata[j,imed] +EGS_Float* EGS_AdvancedApplication::getNbFdata(int j, int imed) { + return the_nist_brems->nb_fdata[imed][j]; +} + +//return the_nist_brems->nb_wdata[j,imed] +EGS_Float* EGS_AdvancedApplication::getNbWdata(int j, int imed) { + return the_nist_brems->nb_wdata[imed][j]; +} + +//return the_nist_brems->nb_idata[j,imed] +int* EGS_AdvancedApplication::getNbIdata(int j, int imed) { + return the_nist_brems->nb_idata[imed][j]; +} + +//return value of the_brempr->delcm[imed] +EGS_Float EGS_AdvancedApplication::getDelcm(int imed) { + return the_brempr->delcm[imed]; +} + +//functions below return value of the_brempr->dl1...dl6[i,imed] +EGS_Float EGS_AdvancedApplication::getDl1(int i, int imed) { + return the_brempr->dl1[imed][i]; +} +EGS_Float EGS_AdvancedApplication::getDl2(int i, int imed) { + return the_brempr->dl2[imed][i]; +} +EGS_Float EGS_AdvancedApplication::getDl3(int i, int imed) { + return the_brempr->dl3[imed][i]; +} +EGS_Float EGS_AdvancedApplication::getDl4(int i, int imed) { + return the_brempr->dl4[imed][i]; +} +EGS_Float EGS_AdvancedApplication::getDl5(int i, int imed) { + return the_brempr->dl5[imed][i]; +} +EGS_Float EGS_AdvancedApplication::getDl6(int i, int imed) { + return the_brempr->dl6[imed][i]; +} +//EGSnrc Mortran calls used in radiative splitting +void EGS_AdvancedApplication::callBrems() { + F77_OBJ_(brems,BREMS)(); +} +void EGS_AdvancedApplication::callAnnih() { + F77_OBJ_(annih,ANNIH)(); +} +void EGS_AdvancedApplication::callAnnihAtRest() { + F77_OBJ_(annih_at_rest,ANNIH_AT_REST)(); +} +void EGS_AdvancedApplication::callPhoto() { + F77_OBJ_(photo,PHOTO)(); +} +void EGS_AdvancedApplication::callPair() { + F77_OBJ_(pair,PAIR)(); +} +void EGS_AdvancedApplication::callCompt() { + F77_OBJ_(compt,COMPT)(); +} +void EGS_AdvancedApplication::callEgsRayleighSampling(int imed, EGS_Float E, EGS_Float gle, EGS_I32 lgle, EGS_Float& costhe, EGS_Float& sinthe) { + int f_imed=imed+1; + F77_OBJ_(egs_rayleigh_sampling,EGS_RAYLEIGH_SAMPLING)(&f_imed,&E,&gle,&lgle,&costhe,&sinthe); +} +EGS_Float EGS_AdvancedApplication::callAliasSample1(int mxbrxs, EGS_Float* nb_xdata, EGS_Float* nb_fdata, EGS_Float* nb_wdata, int* nb_idata) { + return F77_OBJ_(alias_sample1,ALIAS_SAMPLE1)(mxbrxs,nb_xdata,nb_fdata,nb_wdata,nb_idata); +} +//the following definitions are dependent on some array size macros +//having been defined +#ifndef MXSTACK +int EGS_AdvancedApplication::getMxstack() { + return 0; +} +#else +int EGS_AdvancedApplication::getMxstack() { + return MXSTACK; +} +#endif +#ifndef MXBRES +int EGS_AdvancedApplication::getMxbres() { + return 0; +} +#else +int EGS_AdvancedApplication::getMxbres() { + return MXBRES; +} +#endif +#ifndef MXBRXS +int EGS_AdvancedApplication::getMxbrxs() { + return 0; +} +#else +int EGS_AdvancedApplication::getMxbrxs() { + return MXBRXS; +} +#endif + +//get max. energy of source from source +EGS_Float EGS_AdvancedApplication::getEmax() { + return source->getEmax(); +} + +//get value of ZBRANG for medium +EGS_Float EGS_AdvancedApplication::getZbrang(int imed) { + return the_brempr->zbrang[imed]; +} + extern __extc__ void egsHowfar() { CHECK_GET_APPLICATION(app,"egsHowfar()"); int np = the_stack->np-1; @@ -1338,6 +1605,7 @@ extern __extc__ void egsAusgab(EGS_I32 *iarg) { app->top_p.latch = the_stack->latch[np]; app->Np = the_stack->np-1; *iarg = app->userScoring(*iarg); + if (the_stack->wt[np] == 0) return; //allow code to force return to shower } extern __extc__ void egsStartParticle() { @@ -1350,6 +1618,5 @@ extern __extc__ void egsStartParticle() { } the_epcont->idisc = 0; the_useful->medium = app->getMedium(ir)+1; - //egsInformation("start particle: ir=%d medium=%d\n",ir,the_useful->medium); app->startNewParticle(); } diff --git a/HEN_HOUSE/egs++/egs_advanced_application.h b/HEN_HOUSE/egs++/egs_advanced_application.h index b77fc6093..3177f7e57 100644 --- a/HEN_HOUSE/egs++/egs_advanced_application.h +++ b/HEN_HOUSE/egs++/egs_advanced_application.h @@ -227,6 +227,51 @@ class APP_EXPORT EGS_AdvancedApplication : public EGS_Application { int getNp(); int getNpOld(); + // Utility functions for egs_radiative_splitting ausgab objects + //************************************************************ + EGS_Float getRngUniform(); + void getRngAzimuth(EGS_Float &cphi, EGS_Float &sphi); + void addParticleToStack(EGS_Particle p, EGS_Float dnear); + void updateParticleOnStack(int ip, EGS_Particle p, EGS_Float dnear); + EGS_Float getDnear(int np); + int getNmed(); + void setNpold(int npold); + void deleteParticleFromStack(int ip); + EGS_Particle getParticleFromStack(int ip); + EGS_Float getGle(); + int getLgle(EGS_Float gle, int med); + int getIbrdst(); + int getIbcmp(); + EGS_Float getAp(int imed); + int getIbrnist(); + EGS_Float getNbLemin(int imed); + EGS_Float getNbDlei(int imed); + EGS_Float getNbEmin(int imed); + EGS_Float* getNbXdata(int j, int imed); + EGS_Float* getNbFdata(int j, int imed); + EGS_Float* getNbWdata(int j, int imed); + int* getNbIdata(int j, int imed); + EGS_Float getEmax(); + EGS_Float getZbrang(int imed); + EGS_Float getDelcm(int imed); + EGS_Float getDl1(int i, int imed); + EGS_Float getDl2(int i, int imed); + EGS_Float getDl3(int i, int imed); + EGS_Float getDl4(int i, int imed); + EGS_Float getDl5(int i, int imed); + EGS_Float getDl6(int i, int imed); + void callBrems(); + void callAnnih(); + void callAnnihAtRest(); + void callPhoto(); + void callPair(); + void callCompt(); + void callEgsRayleighSampling(int imed, EGS_Float e, EGS_Float gle, EGS_I32 lgle, EGS_Float& costhe, EGS_Float& sinthe); + EGS_Float callAliasSample1(int mxbrxs, EGS_Float* nb_xdata, EGS_Float* nb_fdata, EGS_Float* nb_wdata, int* nb_idata); + int getMxstack(); + int getMxbres(); + int getMxbrxs(); + /* Needed by some sources */ EGS_Float getRM(); /* Turn ON/OFF radiative splitting */ diff --git a/HEN_HOUSE/egs++/egs_application.h b/HEN_HOUSE/egs++/egs_application.h index c7842f38b..6fad85282 100644 --- a/HEN_HOUSE/egs++/egs_application.h +++ b/HEN_HOUSE/egs++/egs_application.h @@ -1190,6 +1190,112 @@ class EGS_EXPORT EGS_Application { //************************************************************ virtual void setLatch(int latch) {}; + //************************************************************ + // Utility functions for egs_radiative_splitting ausgab objects + //************************************************************ + virtual EGS_Float getRngUniform() { + return 0.0; + } + virtual void getRngAzimuth(EGS_Float &cphi, EGS_Float &sphi) {}; + virtual void addParticleToStack(EGS_Particle p, EGS_Float dnear) {}; + virtual EGS_Float getDnear(int np) { + return 0.0; + } + virtual int getNmed() { + return 0; + } + virtual void setNpold(int npold) {}; + virtual void deleteParticleFromStack(int ip) {}; + virtual EGS_Particle getParticleFromStack(int ip) { + EGS_Particle p; + return p; + }; + virtual void updateParticleOnStack(int ip, EGS_Particle p, EGS_Float dnear) {}; + virtual EGS_Float getGle() { + return 0.0; + } + virtual int getLgle(EGS_Float gle, int med) { + return 0; + } + virtual int getIbrdst() { + return 0; + } + virtual int getIbcmp() { + return 0; + } + virtual EGS_Float getAp(int imed) { + return 0.0; + } + virtual int getIbrnist() { + return 0.0; + } + virtual EGS_Float getNbLemin(int imed) { + return 0.0; + } + virtual EGS_Float getNbDlei(int imed) { + return 0.0; + } + virtual EGS_Float getNbEmin(int imed) { + return 0.0; + } + virtual EGS_Float* getNbXdata(int j, int imed) { + return NULL; + } + virtual EGS_Float* getNbFdata(int j, int imed) { + return NULL; + } + virtual EGS_Float* getNbWdata(int j, int imed) { + return NULL; + } + virtual int* getNbIdata(int j, int imed) { + return NULL; + } + virtual EGS_Float getEmax() { + return 0.0; + } + virtual EGS_Float getZbrang(int imed) { + return 0.0; + } + virtual EGS_Float getDelcm(int imed) { + return 0.0; + } + virtual EGS_Float getDl1(int i, int imed) { + return 0.0; + } + virtual EGS_Float getDl2(int i, int imed) { + return 0.0; + } + virtual EGS_Float getDl3(int i, int imed) { + return 0.0; + } + virtual EGS_Float getDl4(int i, int imed) { + return 0.0; + } + virtual EGS_Float getDl5(int i, int imed) { + return 0.0; + } + virtual EGS_Float getDl6(int i, int imed) { + return 0.0; + } + virtual void callBrems() {}; + virtual void callAnnih() {}; + virtual void callAnnihAtRest() {}; + virtual void callPhoto() {}; + virtual void callPair() {}; + virtual void callCompt() {}; + virtual void callEgsRayleighSampling(int imed, EGS_Float e, EGS_Float gle, EGS_I32 lgle, EGS_Float& costhe, EGS_Float& sinthe) {}; + virtual EGS_Float callAliasSample1(int mxbrxs, EGS_Float* nb_xdata, EGS_Float* nb_fdata, EGS_Float* nb_wdata, int* nb_idata) { + return 0.0; + } + virtual int getMxstack() { + return 0; + } + virtual int getMxbres() { + return 0; + } + virtual int getMxbrxs() { + return 0; + } }; #define APP_MAIN(app_name) \ diff --git a/HEN_HOUSE/egs++/egs_ausgab_object.h b/HEN_HOUSE/egs++/egs_ausgab_object.h index ceeda3275..eb43671f3 100644 --- a/HEN_HOUSE/egs++/egs_ausgab_object.h +++ b/HEN_HOUSE/egs++/egs_ausgab_object.h @@ -109,6 +109,9 @@ class EGS_EXPORT EGS_AusgabObject : public EGS_Object { /*! \brief Set the current event */ virtual void setCurrentCase(EGS_I64 ncase) {}; + /*! \brief Initialize variables at the beginning of each shower */ + virtual void initializeData() {}; + /*! \brief Get a short description of this ausgab object. * * Derived classes should set #description to a short diff --git a/HEN_HOUSE/interface/egs_c_interface2.macros b/HEN_HOUSE/interface/egs_c_interface2.macros index 5a1a73151..042c78276 100644 --- a/HEN_HOUSE/interface/egs_c_interface2.macros +++ b/HEN_HOUSE/interface/egs_c_interface2.macros @@ -360,7 +360,11 @@ REPLACE {$CALL-HOWFAR-IN-PHOTON;} WITH {; REPLACE {$AUSCALL(#);} WITH { iarg={P1}; - IF (IAUSFL(IARG+1).NE.0) [call egs_ausgab(iarg);] + IF (IAUSFL(IARG+1).NE.0) [ + call egs_ausgab(iarg); + "allow app to return to shower early" + IF( wt(np) = 0 ) [ np = np-1; return; ] + ] }; diff --git a/HEN_HOUSE/interface/egs_interface2.c b/HEN_HOUSE/interface/egs_interface2.c index 09b5843f3..1c16af8ea 100644 --- a/HEN_HOUSE/interface/egs_interface2.c +++ b/HEN_HOUSE/interface/egs_interface2.c @@ -54,6 +54,9 @@ extern __extc__ struct EGS_IO F77_OBJ_(egs_io,EGS_IO); extern __extc__ struct EGS_VarianceReduction F77_OBJ_(egs_vr,EGS_VR); extern __extc__ struct EGS_Rayleigh F77_OBJ_(rayleigh_inputs,RAYLEIGH_INPUTS); extern __extc__ struct EGS_emfInputs F77_OBJ_(emf_inputs,EMF_INPUTS); +extern __extc__ struct EGS_NistBrems F77_OBJ_(nist_brems,NIST_BREMS); +extern __extc__ struct EGS_Brempr F77_OBJ_(brempr,BREMPR); +extern __extc__ struct EGS_Photin F77_OBJ_(photin,PHOTIN); struct EGS_Stack *the_stack = & F77_OBJ(stack,STACK); struct EGS_Bounds *the_bounds = & F77_OBJ(bounds,BOUNDS); @@ -67,6 +70,9 @@ struct EGS_IO *the_egsio = & F77_OBJ_(egs_io,EGS_IO); struct EGS_VarianceReduction *the_egsvr = & F77_OBJ_(egs_vr,EGS_VR); struct EGS_Rayleigh *the_rayleigh = & F77_OBJ_(rayleigh_inputs,RAYLEIGH_INPUTS); struct EGS_emfInputs *the_emf = & F77_OBJ_(emf_inputs,EMF_INPUTS); +struct EGS_NistBrems *the_nist_brems = & F77_OBJ_(nist_brems,NIST_BREMS); +struct EGS_Brempr *the_brempr = & F77_OBJ_(brempr,BREMPR); +struct EGS_Photin *the_photin = & F77_OBJ_(photin,PHOTIN); extern __extc__ void F77_OBJ_(egs_init_f,EGS_INIT_F)(); extern __extc__ void F77_OBJ(electr,ELECTR)(int *); diff --git a/HEN_HOUSE/interface/egs_interface2.h b/HEN_HOUSE/interface/egs_interface2.h index 6baeca119..a0f844415 100644 --- a/HEN_HOUSE/interface/egs_interface2.h +++ b/HEN_HOUSE/interface/egs_interface2.h @@ -557,7 +557,7 @@ struct EGS_XOptions { /*! \brief A structure corresponding to the \c egs_vr common block The \c egs_vr common block contains variables that can turn on/off - internally implemented variabce reduction techniques. + internally implemented variance reduction techniques. */ struct EGS_VarianceReduction { /*! Maximum energy for which it is allowed to use range rejection. @@ -642,11 +642,61 @@ struct EGS_emfInputs { bool emfield_on; }; + +#define MXBRES 100 +#define MXBRXS 50 +/*! \brief A structure corresponding to the \c NIST_BREMS common block. + + This common block contains variables used to calculate NIST brems cross- + sections + */ +struct EGS_NistBrems { + EGS_Float nb_fdata[MXMED][MXBRES][MXBRXS+1], nb_xdata[MXMED][MXBRES][MXBRXS+1], + nb_wdata[MXMED][MXBRES][MXBRXS]; + EGS_I32 nb_idata[MXMED][MXBRES][MXBRXS]; + EGS_Float nb_emin[MXMED],nb_emax[MXMED],nb_lemin[MXMED],nb_lemax[MXMED], + nb_dle[MXMED],nb_dlei[MXMED],log_ap[MXMED]; +}; + +#define MXEL 50 +#define MXPWR2I 50 +/*! \brief A structure corresponding to the \c BREMPR common block. + + This common block contains variables used to calculate NRC brems cross- + sections + */ +struct EGS_Brempr { + EGS_Float dl1[MXMED][8],dl2[MXMED][8],dl3[MXMED][8],dl4[MXMED][8],dl5[MXMED][8], + dl6[MXMED][8],alphi[MXMED][2],bpar[MXMED][2],delpos[MXMED][2], + wa[MXEL][MXMED],pz[MXEL][MXMED],zelem[MXEL][MXMED],rhoz[MXEL][MXMED], + pwr2i[MXPWR2I],delcm[MXMED],zbrang[MXMED],lzbrang[MXMED]; + EGS_I32 nne[MXMED]; + char asym[2][MXEL][MXMED][4]; +}; + +#define MXRAYFF 100 +#define MXSGE 400 +#define MXGE 2000 +/*! \brief A structure corresponding to the \c PHOTIN common block. + + This common block contains variables used to calculate data used in + photon transport +*/ +struct EGS_Photin { + EGS_Float ebinda[MXMED],ge0[MXMED],ge1[MXMED],gmfp1[MXMED][MXGE],gmfp2[MXMED][MXGE], + gbr10[MXMED][MXGE],gbr11[MXMED][MXGE],gbr20[MXMED][MXGE],gbr21[MXMED][MXGE], + rco0[MXMED],rco1[MXMED],rsct0[MXMED][MXRAYFF],rsct1[MXMED][MXRAYFF], + cohe0[MXMED][MXGE],cohe1[MXMED][MXGE],photonuc0[MXMED][MXGE],photonuc1[MXMED][MXGE], + dpmfp; + EGS_I32 mpgem[MXMED][MXSGE],ngr[MXMED]; +}; + + /*! \brief The address of the mortran \c STACK common block as a pointer to a C-structure of type EGS_Stack */ extern __extc__ struct EGS_Stack *the_stack; -/*! \brief The address of the morrtan \c BOUNDS common block as a +/*! \brief The address of the mortan \c BOUNDS common block as a pointer to a C-structure of type EGS_Bounds */ extern __extc__ struct EGS_Bounds *the_bounds; @@ -690,6 +740,18 @@ extern __extc__ struct EGS_Rayleigh *the_rayleigh; pointer to a C-structure of type EGS_emfInputs */ extern __extc__ struct EGS_emfInputs *the_emf; +/*! \brief The address of the mortran NIST-BREMS common block as a + pointer to a C-structure of type EGS_NistBrems */ +extern __extc__ struct EGS_NistBrems *the_nist_brems; + +/*! \brief The address of the mortran BREMPR common block as a + pointer to a C-structure of type EGS_Brempr */ +extern __extc__ struct EGS_Brempr *the_brempr; + +/*! \brief The address of the mortran PHOTIN common block as a + pointer to a C-structure of type EGS_Photin */ +extern __extc__ struct EGS_Photin *the_photin; + /* ******************* EGSnrc interface functions *************************/