The tutor7pp application is a re-implementation of the original tutor7 mortran application in C++ using EGS_AdvancedApplication. It calculates the deposited, transmitted and reflected energy fractions in addition to the pulse height distribution in requested geometry regions. However, unlike the original tutor7 application, which can only calculate these quantites for a slab geometry for a pencil beam incident at a user defined angle, tutor7pp can perform the simulation in any geometry that can be defined using the egspp geometry package for any particle source that can be defined using the egspp source package. In addition, tutor7pp provides the functionality of being able to restart a calculation, run a calculation in parallel and automatically recombine the results, limit the simulation time to a user defined maximum time, or abort the simulation if a user defined statistical uncertainty has been reached. Despite this added functionality, the tutor7pp source code is much shorter than the original tutor7 implementation in mortran. The reader is encouraged to study the guidelines for developing EGSnrc C++ applications provided in the EGS_Application class documentation before reading this example.
To implement the required functionality we derive an application class from EGS_AdvancedApplication
#include <cstdlib>
using namespace std;
The various include files are needed to gain access to the declarations of the classes and functions that will be used. We then declare various private data members:
int nreg;
int nph;
double Etot;
int rr_flag;
EGS_Float current_weight;
bool deflect_brems;
EGS_Float *ph_de;
int *ph_regions;
static string revision;
Their function is briefly explained by the comments and will become more clear in what follows. We now need a constructor for our application class that takes the command line arguments as input, passes them to the base class constructor and initializes its data
public:
Tutor7_Application(int argc, char **argv) :
nreg(0), nph(0), Etot(0), rr_flag(0), current_weight(1), deflect_brems(false) { };
The call to the
EGS_AdvancedApplication constructor checks the arguments passed on the command line and sets up the name of the application, the input file name, the pegs file name, the output file name, the working directory, batch vs interactive run, parallel job execution (i_parallel and n_parallel) and reads in the input file into an
EGS_Input object (a pointer to this object is available with the protected data member
EGS_Application::input). To demonstrate good coding habits, we declare a destructor for our application that deallocates the memory used
~Tutor7_Application() {
if (score) {
delete score;
}
if (eflu) {
delete eflu;
}
if (gflu) {
delete gflu;
}
if (nph > 0) {
for (int j=0; j<nph; j++) {
delete pheight[j];
}
delete [] pheight;
delete [] ph_regions;
delete [] ph_de;
}
};
We now declare the various virtual functions that we wish to re-implement: the function to describe our user code,
void describeUserCode() const;
the function to initialize the scoring variables,
the function for scoring the quantities of interest at run time,
the functions needed for restarting simulations and parallel runs,
int outputData();
int readData();
void resetCounter();
int addState(istream &data);
the function to output the results of the simulation,
the function that makes the current result known to the
run control object,
void getCurrentResult(double &sum, double &sum2, double &norm,
double &count);
and the protected function called before the simulation of a new particle
protected:
int startNewShower();
};
This completes our application class declaration.
We now provide implementations of the above functions. First we have to declare the static data member revision:
string Tutor7_Application::revision = " ";
Note that the revision string was automatically updated when the code was under CVS revision control. Since the system has been moved to GitHub, however, the revision number is now blank, pending the implementation of a sensible/automatic numbering scheme. The
describeUserCode() function simply outputs the name of our application and the revision numbers of the base application class and our application class using egsInformation:
extern "C" void F77_OBJ_(egs_scale_xcc,EGS_SCALE_XCC)(const int *,const EGS_Float *);
extern "C" void F77_OBJ_(egs_scale_bc,EGS_SCALE_BC)(const int *,const EGS_Float *);
void Tutor7_Application::describeUserCode() const {
"\n ***************************************************"
"\n * *"
"\n * tutor7pp *"
"\n * *"
"\n ***************************************************"
"\n\n");
" EGS_AdvancedApplication %s\n\n",
}
This is useful for determining which version was used to obtain a given result from the logfile produced by the application. The egsSimplifyCVSKey function used above removes the $'s from CVS keys.
The next step is to initialize the scoring variables in initScoring(). We will need a scoring array object to collect energy depositions in all regions of the geometry and the outside region (divided into energy transmitted, defined as the energy of particles exiting the geometry and moving forward, and energy reflected, defined as the energy of particles exiting the geometry and moving backwards) from which we can calculate energy deposition fractions. For this purpose we obtain the number of regions in the geometry
int Tutor7_Application::initScoring() {
nreg = geometry->regions();
and initialize a scoring array with
nreg+2 regions pointed to by our class variable
score: We want the user to be able to specify the geometry regions for which to calculate a pulse height distribution (PHD) by including
:start scoring options:
pulse height regions = 1 5 77
pulse height bins = 125 50 500
:stop scoring options:
in their input file, with the pulse height regions key specifying all regions for which a PHD calculation is requested and the pulse height bins key the number of bins to use for each PHD region (obviously, the number of integer inputs to these two keys must be the same). To implement this functionality, we interogate the input object , which contains all input found in the input file, for a composite key named scoring options
and if such property exists (
i.e.
options is not
null), check for the
pulse height regions and
pulse height bins keys:
vector<string> tmp;
int err = scale->getInput("scale xcc",tmp);
if (!err && tmp.size() == 2) {
if (imed > 0) {
EGS_Float fac = atof(tmp[1].c_str());
egsInformation(
"\n ***** Scaling xcc of medium %d with %g\n",imed,fac);
F77_OBJ_(egs_scale_xcc,EGS_SCALE_XCC)(&imed,&fac);
}
}
delete scale;
}
vector<string> tmp;
int err = scale->getInput("scale bc",tmp);
if (!err && tmp.size() == 2) {
if (imed > 0) {
EGS_Float fac = atof(tmp[1].c_str());
egsInformation(
"\n ***** Scaling bc of medium %d with %g\n",imed,fac);
F77_OBJ_(egs_scale_bc,EGS_SCALE_BC)(&imed,&fac);
}
}
delete scale;
}
vector<string> choices;
choices.push_back("no");
choices.push_back("yes");
deflect_brems = options->
getInput(
"deflect electron after brems",choices,0);
if (deflect_brems) {
egsInformation(
"\n *** Using electron deflection in brems events\n\n");
setAusgabCall(AfterBrems,true);
}
int n_rr;
if (!options->
getInput(
"Russian Roulette",n_rr) && n_rr > 1) {
setAusgabCall(BeforeBrems,true);
setAusgabCall(AfterBrems,true);
setAusgabCall(BeforeAnnihFlight,true);
setAusgabCall(AfterAnnihFlight,true);
setAusgabCall(BeforeAnnihRest,true);
setAusgabCall(AfterAnnihRest,true);
egsInformation(
"\nUsing Russian Roulette with survival probability 1/%d\n",n_rr);
}
vector<int> regions;
int err = options->
getInput(
"pulse height regions",regions);
vector<int> nbins;
int err1 = options->
getInput(
"pulse height bins",nbins);
If the keys exist (
i.e.
err and
err1 are zero), we check that the number of inputs is the same, or, alternatively there is only a single input to the
pulse height bins key (in this case we use the same number of bins for all PHDs), and if this is not satisfied issue a warning so that the user is made aware of a potential mistake in the input file,
if (!err && !err1) {
if (regions.size() != nbins.size() && nbins.size() != 1)
egsWarning(
"initScoring(): you must input the same "
"number of 'regions' and 'bins' inputs or a single 'bins'"
" input\n");
otherwise we have to initialize scoring arrays for the requested PHDs. We first alocate a temporary array of pointers to scoring objects and set all of them to
null: else {
for (int i=0; i<nreg+2; i++) {
tmp[i] = 0;
We then loop over the number of PHD regions
}
for (int j=0; j<regions.size(); j++) {
For each requested PHD region we set the number of bins (possibly warning the user, if the number is less than 1),
int nb = nbins.size() == 1 ? nbins[0] : nbins[j];
if (nb < 1) {
egsWarning(
"zero bins for region %d?\n",regions[j]);
and check that the region index is valid
}
if (regions[j] < -1 || regions[j] > nreg) {
egsWarning(
"invalid region index %d\n",regions[j]);
If the number of bins is positive and the region index is valid, we check if there is not already a PHD scoring object for this region initilized and warn the user if there is
}
if (nb > 0 && regions[j] >= 0 && regions[j] < nreg+2) {
int ij = regions[j];
"PHD object in region %d => ignoring it\n",ij);
(perhaps the user has made a mistake and input the same PHD region twice), otherwise we construct a scoring array with
nb bins, increment the number of requested PHDs by one and store the pointer to the newly constructed scoring array in
tmp: When the loop over requested PHD regions is completed,
nph has the number of valid PHD requests and the elements of
tmp that are not
null have the pointers to the scorring array objects allocated for the calculation of the PHDs. We then allocate an array of pointers to scorring array objects with the required size,
an array with region indeces for each requested PHD,
ph_regions = new int [nph];
and an array with bin widths for each PHD
ph_de = new EGS_Float [nph];
We then obtain the maximum source energy
EGS_Float Emax = source->getEmax();
and loop over all regions
int iph = 0;
for (int j=0; j<nreg+2; j++) {
If
tmp[j] is not
null, we set the
iph't element of the just allocated arrays to contain the pointer to the scorring array, the region index of the PHD, and the bin width:
if (tmp[j]) {
pheight[iph] = tmp[j];
ph_regions[iph] = j;
int nbin = pheight[iph]->
bins();
ph_de[iph++] = Emax/nbin;
}
}
}
We then delete
tmp to avoid memory leaks
This completes the initilization of PHD scoring. If the user has provided
scoring options input but the input is incorrect (indicated by
err and/or
err1 not being zero), we warn them so that they can check their input file
}
else egsWarning(
"initScoring(): you must provide both, 'regions'"
" and 'bins' input\n");
We then delete the
scoring options input as it is no longer needed
delete options;
}
return 0;
}
(Note: if we were expecting that someone might want to extend the functionality of our tutor7pp application by deriving from the Tutor7_Application class, we should have introduced a protected data member to point to the
scroring options property and should not delete it here. In this case we also should have split the class declaration into a separate header file). This completes the
initScoring() implementation.
The implementation of our ausgab() function is very simple: for enertgy deposition events (i.e. iarg<=4) we determine the region index of the top particle on the particle stack and score the energy deposited weighted with the statistical weight of the particle into the corresponding scoring region using the score() method of our scoring array object:
int Tutor7_Application::ausgab(int iarg) {
if (iarg <= 4) {
ir = nreg+1;
}
if (aux > 0) {
score->score(ir,aux);
}
The calculation of requested pulse height distributions will be done in the
startNewShower() function (see below).
The outputData() function must store the accumulated results to a data stream (click on outputData() below to see the base class documentation of this virtual function). We use the base class implementation to store the data of the base application and return the generated error code, if an error occured:
if (ir == nreg+1) {
int bin = (int)(sqrt(r2)*10.);
if (bin < 200) {
if (aux > 0) {
}
}
}
return 0;
}
if (iarg == BeforeBrems || iarg == BeforeAnnihRest || (iarg == BeforeAnnihFlight &&
rr_flag = 1;
return 0;
}
if (iarg == AfterBrems && deflect_brems) {
EGS_Float beta = sqrt(tau*(tau+2))/(tau+1);
EGS_Float eta = 2*rndm->getUniform()-1;
EGS_Float cost = (beta + eta)/(1 + beta*eta);
EGS_Float sint = 1 - cost*cost;
if (sint > 0) {
sint = sqrt(sint);
EGS_Float cphi, sphi;
rndm->getAzimuth(cphi,sphi);
u.rotate(cost,sint,cphi,sphi);
}
}
if (iarg == AfterBrems || iarg == AfterAnnihRest || iarg == AfterAnnihFlight) {
if (iarg == AfterBrems && rr_flag) {
}
rr_flag = 0;
return 0;
}
return 0;
}
int Tutor7_Application::outputData() {
if (err) {
return err;
We then write our own data to the data stream pointed to by
data_out: the total energy deposited so far
(*data_out) << " " << Etot << endl;
the energy fractions accumulated in the scorring array
if (!score->storeState(*data_out)) {
and the pulse height distribution results
return 101;
}
for (int j=0; j<nph; j++) {
if (!pheight[j]->storeState(*data_out)) {
return 102+j;
}
Note that
data_out is set up to point to the output stream created by opening the
.egsdat file for writing in the base application class.
The function directly innverse to outputData() is readData(). The implementation is therefore basically the same as outputData() except that now we read the data from the stream pointed to by data_in:
}
if (!eflu->storeState(*data_out)) {
return 301;
}
if (!gflu->storeState(*data_out)) {
return 302;
}
return 0;
}
int Tutor7_Application::readData() {
if (err) {
return err;
(*data_in) >> Etot;
if (!score->setState(*data_in)) {
return 101;
}
for (int j=0; j<nph; j++) {
if (!pheight[j]->setState(*data_in)) {
return 102+j;
}
The resetCounter() function, which is needed for combining parallel runs and is called before the loop over available data files, must put our application class into a state with zero particles simulated. We therefore use the base class function to reset the data of the base application class,
}
if (!eflu->setState(*data_in)) {
return 301;
}
if (!gflu->setState(*data_in)) {
return 302;
}
return 0;
}
void Tutor7_Application::resetCounter() {
set the total energy deposited to zero and use the
reset() functions of the scoring arrays to reset these:
Etot = 0;
for (int j=0; j<nph; j++) {
pheight[j]->reset();
}
The addState() function is used to combine the results of parallel runs. It is therefore similar to the readData() function except that now the data read from the input stream is added to the existing data instead of replacing the state of the application with these data. We therefore use the base class addState() function,
eflu->reset();
gflu->reset();
}
int Tutor7_Application::addState(istream &data) {
if (err) {
return err;
add the total energy deposited into the geometry,
construct a temporary scoring array object with
nreg+2 regions,
data >> etot_tmp;
Etot += etot_tmp;
read the data into this object,
and add it to the existing data
return 101;
}
(*score) += tmp;
We do the same for the pulse height distributions:
for (int j=0; j<nph; j++) {
return 102 + j;
}
(*pheight[j]) += tmpj;
}
The outputResults() function must output the results of a simulation into the logfile (or the standard output, if the application is run interactively). In our implementation we first inform the user about the last statistically independent history simulated and the total energy deposited into the geometry:
if (!tmp1.setState(data)) {
return 301;
}
(*eflu) += tmp1;
if (!tmp1.setState(data)) {
return 302;
}
(*gflu) += tmp1;
return 0;
}
void Tutor7_Application::outputResults() {
(int)current_case,Etot);
We then use the
reportResults() method of the energy fractions scoring array to report the deposited energy fractions:
double norm = ((double)current_case)/Etot;
egsInformation(
"\n\n======================================================\n");
Because the scoring array object divides the results by the number of statistically independent showers (
last_case) but we want the normalization to be fraction of the energy deposited, we set the normalization constant to be
last_case/Etot. The string passed as a second argument will be output as a title and the
false argument indicates that we want the uncertainties to be absolute instead of relative as percent of the local result. The last argument is a format string used to output the result. If omitted or set to
null, a default format string will be used (using a g format for outputing the floating numbers). If the user has requested the calculation of one or more pulse height distributions (
i.e.
nph>0), we output a title depending on the number of distributions calculated:
egsInformation(
"======================================================\n");
and then in a loop over the number of pulse height distributions calculated
egsInformation(
"The first and last items in the following list of energy fractions are the reflected and transmitted energy, respectively. These two values are only meaningful if the source is directed in the positive z-direction. The remaining values are the deposited energy fractions in the regions of the geometry, but notice that the identifying index is the region number offset by 1 (ir+1).");
we have a loop over the number of bins for the pulse height distribution in this region, obtain the result from the scoring array and output a data triplet consisting of the midpoint energy of the bin, the result in the bin and the uncertainty of the result:
score->reportResults(norm,
"ir+1 | Reflected, deposited, or transmitted energy fraction",false,
" %d %12.6e +/- %12.6e %c\n");
if (nph > 0) {
if (nph > 1) {
egsInformation(
"\n\n======================================================\n");
"======================================================\n\n");
}
else {
"======================================================\n\n",
ph_regions[0]);
}
for (int j=0; j<nph; j++) {
"----------------\n\n",ph_regions[j]);
double f,df;
for (
int i=0; i<pheight[j]->bins(); i++) {
pheight[j]->currentResult(i,f,df);
f/ph_de[j],df/ph_de[j]);
}
}
Because the scoring array object already divides by the number of statistically independent events, we only divide by the bin width to have the results normalized as counts per incident particle per MeV.
The getCurrentResult() function is used by the run control object to provide a single result for our simulation that is the combined result from all jobs in parallel runs. We arbitrarily decide to define the single result of our simulation monitored at run time as the energy fraction reflected (stored in the zeroth region of the score scoring array). Hence, the implementation of this function looks like this:
In the startNewShower() function we must collect the energy being imparted into the geometry and inform our scoring objects when the simulation of a new statistically independent event begins. In addition, we implement the scoring of the requested pulse height distributions (PHD) here. We start by collecting the energy entering the geometry
}
void Tutor7_Application::getCurrentResult(double &sum, double &sum2,
double &norm, double &count) {
count = current_case;
norm = Etot > 0 ? count/Etot : 0;
score->currentScore(0,sum,sum2);
}
int Tutor7_Application::startNewShower() {
Etot += p.E*p.wt;
The parameters of the particle that is about to be simulated are available in the
EGS_Application::p protected data member. We then check if the particle about to be simulated is statistically independent from the last particle simulated by comparing the inherited protected data members
EGS_Application::current_case and
EGS_Application::last_case:
If they are the same, then the two particles belong to the same statistically independent case and we don't need to do anything. This can happen
e.g. when using a
phase space file or a
BEAM simulation source , see
e.g. [Walters
et.al., Med.Phys. 29 (2002) 2745] for a detailed discussion of the correct implementation of a history-by-history statistical analysis. If the user has requested the calculation of PHDs (
i.e.
nph>0), we loop over all requested PHDs
return res;
}
if (current_case != last_case) {
if (nph > 0) {
for (int j=0; j<nph; j++) {
inform the scoring array used to collect this PHD that a new statistically independent event is about to begin,
pheight[j]->setHistory(current_case);
obtain the energy deposited in the PHD region in the last statistically independent shower,
int ireg = ph_regions[j];
EGS_Float edep = score->currentScore(ireg+1);
and, if this energy is greater than zero, determine the PHD bin to which such an energy deposition belongs and add a count to this PHD bin
if (edep > 0) {
int ibin = min((int)(edep/(current_weight*ph_de[j])), pheight[j]->bins()-1);
if (ibin >= 0 && ibin < pheight[j]->bins()) {
pheight[j]->score(ibin,1);
}
}
}
Note that in the above we divide by
current_weight. This is needed because all energy depositions are weighted with the particle weights that may not be constant and also may be different from unity (
e.g. for a
collimated point source ). We remember the statistical weight of the particles intiating the showers in
current_weight (see below). We then inform the energy deposition scoring array that a new statistically independent event begins
}
score->setHistory(current_case);
and set the last simulated case to the current case
eflu->setHistory(current_case);
gflu->setHistory(current_case);
last_case = current_case;
}
and remember the statistical weight of the incident particle.
current_weight = p.wt;
return 0;
}
This completes the implementation of our Tutor7_Application class. One needs more code compared to an application derived from EGS_SimpleApplication (see e.g. this example). However, a closer examination shows that a significant part of the code is neede to implement the restart/parallel runs/combine results functionality and to provide the user with the capability to specify the regions for PHD calculation in the input file. Without this added functionality, the tutor7pp source code would be not much longer than tutor2pp.
We now need a short main function. We construct an instance of the Tutor7_Application class passing the command line arguments to its constructor
#ifdef BUILD_APP_LIB
APP_LIB(Tutor7_Application);
We then initialize the application
and return the error code if some error occured (
e.g. the cross section data for the media in the geometry was not found in the specified PEGS data file, there was a mistake in the specification of the geometry, etc.)
We then run the simulation
and return the error code if some error occured,
otherwise we finish the simulation
That's all.
The final step is to provide a Makefile for our application. We will use the standard makefile for building C++ applications provided with the EGSnrc distribution. To do so, we include the make config files,
include $(EGS_CONFIG)
include $(SPEC_DIR)egspp.spec
include $(SPEC_DIR)egspp_$(my_machine).conf
set the application name
specify that no additional macros are needed to build the EGSnrc mortran back-end,
specify that we are deriving from
EGS_AdvancedApplication,
specify the set of mortran sources to be the predefined set of sources for advanced applications,
the additional header files that our application depends upon,
and finally include the standard set of rules for building C++ applications
Here is the complete source code of the tutor7pp application:
#include <cstdlib>
using namespace std;
int nreg;
int nph;
double Etot;
int rr_flag;
EGS_Float current_weight;
bool deflect_brems;
EGS_Float *ph_de;
int *ph_regions;
static string revision;
public:
Tutor7_Application(int argc, char **argv) :
nreg(0), nph(0), Etot(0), rr_flag(0), current_weight(1), deflect_brems(false) { };
~Tutor7_Application() {
if (score) {
delete score;
}
if (eflu) {
delete eflu;
}
if (gflu) {
delete gflu;
}
if (nph > 0) {
for (int j=0; j<nph; j++) {
delete pheight[j];
}
delete [] pheight;
delete [] ph_regions;
delete [] ph_de;
}
};
void describeUserCode() const;
int initScoring();
int ausgab(int iarg);
int outputData();
int readData();
void resetCounter();
int addState(istream &data);
void outputResults();
void getCurrentResult(double &sum, double &sum2, double &norm,
double &count);
protected:
int startNewShower();
};
string Tutor7_Application::revision = " ";
extern "C" void F77_OBJ_(egs_scale_xcc,EGS_SCALE_XCC)(const int *,const EGS_Float *);
extern "C" void F77_OBJ_(egs_scale_bc,EGS_SCALE_BC)(const int *,const EGS_Float *);
void Tutor7_Application::describeUserCode() const {
"\n ***************************************************"
"\n * *"
"\n * tutor7pp *"
"\n * *"
"\n ***************************************************"
"\n\n");
" EGS_AdvancedApplication %s\n\n",
}
int Tutor7_Application::initScoring() {
nreg = geometry->regions();
if (options) {
vector<string> tmp;
int err = scale->
getInput(
"scale xcc",tmp);
if (!err && tmp.size() == 2) {
if (imed > 0) {
EGS_Float fac = atof(tmp[1].c_str());
egsInformation(
"\n ***** Scaling xcc of medium %d with %g\n",imed,fac);
F77_OBJ_(egs_scale_xcc,EGS_SCALE_XCC)(&imed,&fac);
}
}
delete scale;
}
vector<string> tmp;
int err = scale->
getInput(
"scale bc",tmp);
if (!err && tmp.size() == 2) {
if (imed > 0) {
EGS_Float fac = atof(tmp[1].c_str());
egsInformation(
"\n ***** Scaling bc of medium %d with %g\n",imed,fac);
F77_OBJ_(egs_scale_bc,EGS_SCALE_BC)(&imed,&fac);
}
}
delete scale;
}
vector<string> choices;
choices.push_back("no");
choices.push_back("yes");
deflect_brems = options->
getInput(
"deflect electron after brems",choices,0);
if (deflect_brems) {
egsInformation(
"\n *** Using electron deflection in brems events\n\n");
setAusgabCall(AfterBrems,true);
}
int n_rr;
if (!options->
getInput(
"Russian Roulette",n_rr) && n_rr > 1) {
setAusgabCall(BeforeBrems,true);
setAusgabCall(AfterBrems,true);
setAusgabCall(BeforeAnnihFlight,true);
setAusgabCall(AfterAnnihFlight,true);
setAusgabCall(BeforeAnnihRest,true);
setAusgabCall(AfterAnnihRest,true);
egsInformation(
"\nUsing Russian Roulette with survival probability 1/%d\n",n_rr);
}
vector<int> regions;
int err = options->
getInput(
"pulse height regions",regions);
vector<int> nbins;
int err1 = options->
getInput(
"pulse height bins",nbins);
if (!err && !err1) {
if (regions.size() != nbins.size() && nbins.size() != 1)
egsWarning(
"initScoring(): you must input the same "
"number of 'regions' and 'bins' inputs or a single 'bins'"
" input\n");
else {
for (int i=0; i<nreg+2; i++) {
tmp[i] = 0;
}
for (int j=0; j<regions.size(); j++) {
int nb = nbins.size() == 1 ? nbins[0] : nbins[j];
if (nb < 1) {
egsWarning(
"zero bins for region %d?\n",regions[j]);
}
if (regions[j] < -1 || regions[j] > nreg) {
egsWarning(
"invalid region index %d\n",regions[j]);
}
if (nb > 0 && regions[j] >= 0 && regions[j] < nreg+2) {
int ij = regions[j];
"PHD object in region %d => ignoring it\n",ij);
else {
++nph;
}
}
}
if (nph > 0) {
ph_regions = new int [nph];
ph_de = new EGS_Float [nph];
EGS_Float Emax = source->getEmax();
int iph = 0;
for (int j=0; j<nreg+2; j++) {
if (tmp[j]) {
pheight[iph] = tmp[j];
ph_regions[iph] = j;
int nbin = pheight[iph]->
bins();
ph_de[iph++] = Emax/nbin;
}
}
}
delete [] tmp;
}
}
else egsWarning(
"initScoring(): you must provide both, 'regions'"
" and 'bins' input\n");
delete options;
}
return 0;
}
int Tutor7_Application::ausgab(int iarg) {
if (iarg <= 4) {
ir = nreg+1;
}
if (aux > 0) {
score->score(ir,aux);
}
if (ir == nreg+1) {
int bin = (int)(sqrt(r2)*10.);
if (bin < 200) {
if (aux > 0) {
}
}
}
return 0;
}
if (iarg == BeforeBrems || iarg == BeforeAnnihRest || (iarg == BeforeAnnihFlight &&
rr_flag = 1;
return 0;
}
if (iarg == AfterBrems && deflect_brems) {
EGS_Float beta = sqrt(tau*(tau+2))/(tau+1);
EGS_Float eta = 2*rndm->getUniform()-1;
EGS_Float cost = (beta + eta)/(1 + beta*eta);
EGS_Float sint = 1 - cost*cost;
if (sint > 0) {
sint = sqrt(sint);
EGS_Float cphi, sphi;
rndm->getAzimuth(cphi,sphi);
u.rotate(cost,sint,cphi,sphi);
}
}
if (iarg == AfterBrems || iarg == AfterAnnihRest || iarg == AfterAnnihFlight) {
if (iarg == AfterBrems && rr_flag) {
}
rr_flag = 0;
return 0;
}
return 0;
}
int Tutor7_Application::outputData() {
if (err) {
return err;
}
(*data_out) << " " << Etot << endl;
if (!score->storeState(*data_out)) {
return 101;
}
for (int j=0; j<nph; j++) {
if (!pheight[j]->storeState(*data_out)) {
return 102+j;
}
}
if (!eflu->storeState(*data_out)) {
return 301;
}
if (!gflu->storeState(*data_out)) {
return 302;
}
return 0;
}
int Tutor7_Application::readData() {
if (err) {
return err;
}
(*data_in) >> Etot;
if (!score->setState(*data_in)) {
return 101;
}
for (int j=0; j<nph; j++) {
if (!pheight[j]->setState(*data_in)) {
return 102+j;
}
}
if (!eflu->setState(*data_in)) {
return 301;
}
if (!gflu->setState(*data_in)) {
return 302;
}
return 0;
}
void Tutor7_Application::resetCounter() {
score->reset();
Etot = 0;
for (int j=0; j<nph; j++) {
pheight[j]->reset();
}
eflu->reset();
gflu->reset();
}
int Tutor7_Application::addState(istream &data) {
if (err) {
return err;
}
double etot_tmp;
data >> etot_tmp;
Etot += etot_tmp;
return 101;
}
(*score) += tmp;
for (int j=0; j<nph; j++) {
if (!tmpj.setState(data)) {
return 102 + j;
}
(*pheight[j]) += tmpj;
}
if (!tmp1.setState(data)) {
return 301;
}
(*eflu) += tmp1;
if (!tmp1.setState(data)) {
return 302;
}
(*gflu) += tmp1;
return 0;
}
void Tutor7_Application::outputResults() {
(int)current_case,Etot);
double norm = ((double)current_case)/Etot;
egsInformation(
"\n\n======================================================\n");
egsInformation(
"======================================================\n");
egsInformation(
"The first and last items in the following list of energy fractions are the reflected and transmitted energy, respectively. These two values are only meaningful if the source is directed in the positive z-direction. The remaining values are the deposited energy fractions in the regions of the geometry, but notice that the identifying index is the region number offset by 1 (ir+1).");
score->reportResults(norm,
"ir+1 | Reflected, deposited, or transmitted energy fraction",false,
" %d %12.6e +/- %12.6e %c\n");
if (nph > 0) {
if (nph > 1) {
egsInformation(
"\n\n======================================================\n");
"======================================================\n\n");
}
else {
"======================================================\n\n",
ph_regions[0]);
}
for (int j=0; j<nph; j++) {
"----------------\n\n",ph_regions[j]);
double f,df;
for (
int i=0; i<pheight[j]->bins(); i++) {
pheight[j]->currentResult(i,f,df);
f/ph_de[j],df/ph_de[j]);
}
}
}
}
void Tutor7_Application::getCurrentResult(double &sum, double &sum2,
double &norm, double &count) {
count = current_case;
norm = Etot > 0 ? count/Etot : 0;
score->currentScore(0,sum,sum2);
}
int Tutor7_Application::startNewShower() {
Etot += p.E*p.wt;
if (res) {
return res;
}
if (current_case != last_case) {
if (nph > 0) {
for (int j=0; j<nph; j++) {
pheight[j]->setHistory(current_case);
int ireg = ph_regions[j];
EGS_Float edep = score->currentScore(ireg+1);
if (edep > 0) {
int ibin = min((int)(edep/(current_weight*ph_de[j])), pheight[j]->bins()-1);
if (ibin >= 0 && ibin < pheight[j]->bins()) {
pheight[j]->score(ibin,1);
}
}
}
}
score->setHistory(current_case);
eflu->setHistory(current_case);
gflu->setHistory(current_case);
last_case = current_case;
}
current_weight = p.wt;
return 0;
}
#ifdef BUILD_APP_LIB
APP_LIB(Tutor7_Application);
#else
APP_MAIN(Tutor7_Application);
#endif
