diff --git a/HEN_HOUSE/egs++/Makefile b/HEN_HOUSE/egs++/Makefile index 619662264..6e9851051 100644 --- a/HEN_HOUSE/egs++/Makefile +++ b/HEN_HOUSE/egs++/Makefile @@ -27,6 +27,7 @@ # Frederic Tessier # Manuel Stoeckl # Max Orok +# Alexandre Demelo # ############################################################################### @@ -53,7 +54,8 @@ geometry_libs = egs_planes egs_cd_geometry egs_gtransformed egs_nd_geometry \ egs_cones egs_gstack egs_prism egs_union egs_pyramid egs_conez\ egs_space egs_elliptic_cylinders egs_smart_envelope \ egs_vhp_geometry egs_octree egs_roundrect_cylinders \ - egs_rz egs_autoenvelope egs_lattice egs_glib egs_mesh + egs_rz egs_autoenvelope egs_lattice egs_glib egs_mesh \ + egs_dynamic_geometry source_libs = egs_collimated_source egs_isotropic_source egs_parallel_beam \ egs_point_source egs_source_collection egs_transformed_source \ @@ -64,7 +66,7 @@ source_libs = egs_collimated_source egs_isotropic_source egs_parallel_beam \ shape_libs = egs_circle egs_ellipse egs_extended_shape egs_gaussian_shape \ egs_line_shape egs_polygon_shape egs_rectangle egs_shape_collection \ egs_voxelized_shape egs_spherical_shell egs_conical_shell \ - egs_circle_perpendicular + egs_circle_perpendicular egs_dynamic_shape aobject_libs = egs_track_scoring egs_dose_scoring egs_radiative_splitting \ egs_phsp_scoring egs_fluence_scoring diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.cpp b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.cpp index 958a8b71c..883f6e7b1 100644 --- a/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.cpp +++ b/HEN_HOUSE/egs++/ausgab_objects/egs_phsp_scoring/egs_phsp_scoring.cpp @@ -242,7 +242,7 @@ void EGS_PhspScoring::reportResults() { void EGS_PhspScoring::storeParticle(EGS_I64 ncase) { //if user requested mu scoring, check if mu is available - if (score_mu && app->getMU() < 0) { + if (score_mu && app->getTimeIndex() < 0) { egsWarning("\nEGS_PhspScoring: User requested mu scoring, but mu is inavailable with this source.\n"); egsWarning("Turning off mu scoring.\n"); score_mu=false; @@ -278,7 +278,7 @@ void EGS_PhspScoring::storeParticle(EGS_I64 ncase) { p_stack[phsp_index].w = app->top_p.u.z; p_stack[phsp_index].q = app->top_p.q; if (score_mu) { - p_stack[phsp_index].mu = app->getMU(); + p_stack[phsp_index].mu = app->getTimeIndex(); } p_stack[phsp_index++].latch = app->top_p.latch; diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_track_scoring/egs_track_scoring.cpp b/HEN_HOUSE/egs++/ausgab_objects/egs_track_scoring/egs_track_scoring.cpp index bfcb82721..e6a9c2bb5 100644 --- a/HEN_HOUSE/egs++/ausgab_objects/egs_track_scoring/egs_track_scoring.cpp +++ b/HEN_HOUSE/egs++/ausgab_objects/egs_track_scoring/egs_track_scoring.cpp @@ -24,6 +24,7 @@ # Author: Iwan Kawrakow, 2009 # # Contributors: Georgi Gerganov +# Alexandre Demelo # ############################################################################### */ @@ -41,7 +42,7 @@ EGS_TrackScoring::EGS_TrackScoring(const string &Name, EGS_ObjectFactory *f) : EGS_AusgabObject(Name,f), m_pts(0), m_start(0), m_stop(1024), m_lastCase(-1), m_nScore(0), m_bufSize(16), m_score(false), m_didScore(false), - m_score_photons(true), m_score_electrons(true), m_score_positrons(true), m_fnExtra("") { + m_score_photons(true), m_score_electrons(true), m_score_positrons(true), m_fnExtra(""), m_include_time(false) { otype = "EGS_TrackScoring"; } @@ -77,8 +78,17 @@ void EGS_TrackScoring::setApplication(EGS_Application *App) { sprintf(buf,"_w%d",i_parallel); fname += buf; } - fname += ".ptracks"; - m_pts = new EGS_ParticleTrackContainer(fname.c_str(),m_bufSize); + // if incltime is false use .ptracks. If incltime is true use the new file + // format .syncptracks + if (m_include_time) { + fname += ".syncptracks"; + } + else { + fname += ".ptracks"; + } + // create new particleTrackContainer using the m_include_time boolean which + // controls time index writting and filetype + m_pts = new EGS_ParticleTrackContainer(fname.c_str(),m_bufSize,m_include_time); description = "\nParticle Track Scoring ("; description += name; @@ -90,6 +100,8 @@ void EGS_TrackScoring::setApplication(EGS_Application *App) { description += m_score_electrons ? "YES\n" : "NO\n"; description += " - Scoring positron tracks = "; description += m_score_positrons ? "YES\n" : "NO\n"; + description += " - Include time index = "; + description += m_include_time ? "YES\n" : "NO\n"; description += " - First event to score = "; char buf[32]; sprintf(buf,"%lld\n",m_start); @@ -114,7 +126,6 @@ void EGS_TrackScoring::reportResults() { } } - extern "C" { EGS_TRACK_SCORING_EXPORT EGS_AusgabObject *createAusgabObject(EGS_Input *input, @@ -130,6 +141,11 @@ extern "C" { bool scph = input->getInput("score photons",sc_options,true); bool scel = input->getInput("score electrons",sc_options,true); bool scpo = input->getInput("score positrons",sc_options,true); + + // include time index lets the program know whether to write the time + // index to the tracks file and determines the filetype (ptracks or + // syncptracks); false in argument makes time inclusion false by default + bool incltime = input->getInput("include time index",sc_options,false); if (!scph && !scel && !scpo) { return 0; } @@ -144,6 +160,7 @@ extern "C" { result->setScorePhotons(scph); result->setScoreElectrons(scel); result->setScorePositrons(scpo); + result->setIncludeTime(incltime); // incltime boolean is set from aquired input for the trackscoring object (sets m_include_time) result->setFirstEvent(first); result->setLastEvent(last); result->setBufferSize(bufSize); diff --git a/HEN_HOUSE/egs++/ausgab_objects/egs_track_scoring/egs_track_scoring.h b/HEN_HOUSE/egs++/ausgab_objects/egs_track_scoring/egs_track_scoring.h index d8e3c2a77..f445f5b8f 100644 --- a/HEN_HOUSE/egs++/ausgab_objects/egs_track_scoring/egs_track_scoring.h +++ b/HEN_HOUSE/egs++/ausgab_objects/egs_track_scoring/egs_track_scoring.h @@ -24,6 +24,7 @@ # Author: Iwan Kawrakow, 2009 # # Contributors: Georgi Gerganov +# Alexandre Demelo # ############################################################################### */ @@ -78,6 +79,7 @@ This ausgab object is specified via score photons = yes or no # optional, yes assumed if missing score electrons = yes or no # optional, yes assumed if missing score positrons = yes or no # optional, yes assumed if missing + include time index = yes or no #optional, no assumed if missing start scoring = event_number # optional, 0 assumed if missing stop scoring = event_number # optional, 1024 assumed if missing buffer size = size # optional, 1024 assumed if missing @@ -116,11 +118,13 @@ class EGS_TRACK_SCORING_EXPORT EGS_TrackScoring : public EGS_AusgabObject { } else if (iarg == EGS_Application::BeforeTransport) { int q = app->top_p.q; + EGS_Float timeindex = app->getTimeIndex(); // get current time index from the source (through the application.cpp) if ((q == 0 && m_score_photons) || (q == -1 && m_score_electrons) || (q == 1 && m_score_positrons)) { m_pts->startNewTrack(np); m_pts->setCurrentParticleInfo(new EGS_ParticleTrack::ParticleInfo(q)); + m_pts->setCurrentTimeIndex(timeindex); // set current time index for the particle track container object m_pts->addVertex(np,new EGS_ParticleTrack::Vertex(app->top_p.x,app->top_p.E)); } } @@ -158,6 +162,9 @@ class EGS_TRACK_SCORING_EXPORT EGS_TrackScoring : public EGS_AusgabObject { void setScorePositrons(bool score) { m_score_positrons = score; }; + void setIncludeTime(bool score) { + m_include_time = score; + }; void setFirstEvent(EGS_I64 first) { m_start = first; }; @@ -188,6 +195,7 @@ class EGS_TRACK_SCORING_EXPORT EGS_TrackScoring : public EGS_AusgabObject { bool m_score_photons; //!< Score photon tracks? bool m_score_electrons; //!< Score electron tracks? bool m_score_positrons; //!< Score positron tracks? + bool m_include_time; //!< include time index in tracks file? string m_fnExtra; //!< String to append to output file name diff --git a/HEN_HOUSE/egs++/egs_application.cpp b/HEN_HOUSE/egs++/egs_application.cpp index e638e843f..458e85987 100644 --- a/HEN_HOUSE/egs++/egs_application.cpp +++ b/HEN_HOUSE/egs++/egs_application.cpp @@ -28,6 +28,7 @@ # Blake Walters # Reid Townson # Hubert Ho +# Alexandre Demelo # ############################################################################### */ @@ -47,6 +48,7 @@ #include "egs_base_source.h" #include "egs_simple_container.h" #include "egs_ausgab_object.h" +#include "egs_base_geometry.h" #include #include @@ -925,8 +927,14 @@ int EGS_Application::simulateSingleShower() { " attempts\n"); return 1; } + setTimeIndex(-1); current_case = source->getNextParticle(rndm,p.q,p.latch,p.E,p.wt,p.x,p.u); + + // For dynamic geometries, update positions according to the current + // time index, which may have been set by getNextParticle + geometry->getNextGeom(rndm); + ireg = geometry->isWhere(p.x); if (ireg < 0) { EGS_Float t = veryFar; diff --git a/HEN_HOUSE/egs++/egs_application.h b/HEN_HOUSE/egs++/egs_application.h index c7842f38b..c8bf3a403 100644 --- a/HEN_HOUSE/egs++/egs_application.h +++ b/HEN_HOUSE/egs++/egs_application.h @@ -27,6 +27,7 @@ # Ernesto Mainegra-Hing # Blake Walters # Reid Townson +# Alexandre Demelo # ############################################################################### */ @@ -694,16 +695,23 @@ class EGS_EXPORT EGS_Application { geometry->getLabelRegions(str, regs); } - /*! \brief Returns the value of the \a mu synchronization parameter + /*! \brief Returns the value of the \a time synchronization parameter - The parameter, \a mu, is a random number on \a [0,1) associated with each + The parameter, \a time, is a random number on \a [0,1) associated with each primary history and is retrieved from \a source. It can be used to - synchronize geometric parameters throughout a simulation. If \a mu is - not available in \a source (i.e., the \a getMu function has not been + synchronize geometric parameters throughout a simulation. If \a time is + not available in \a source (i.e., the \a getTime function has not been reimplemented in \a source), then this returns -1. */ - EGS_Float getMU() { - return source->getMu(); + EGS_Float getTimeIndex() { + return source->getTimeIndex(); + } + + /*! Sets the value of the time synchronization parameter. This will mainly + * be used by the dynamic geometry if it does not receive time from a source + * */ + void setTimeIndex(EGS_Float temp_time) { + source->setTimeIndex(temp_time); } /*! \brief User scoring function for accumulation of results and VRT implementation diff --git a/HEN_HOUSE/egs++/egs_base_geometry.h b/HEN_HOUSE/egs++/egs_base_geometry.h index 9bfed973b..141409e9e 100644 --- a/HEN_HOUSE/egs++/egs_base_geometry.h +++ b/HEN_HOUSE/egs++/egs_base_geometry.h @@ -29,6 +29,7 @@ # Reid Townson # Ernesto Mainegra-Hing # Hugo Bouchard +# Alexandre Demelo # ############################################################################### */ @@ -43,6 +44,7 @@ #define EGS_BASE_GEOMETRY_ #include "egs_vector.h" +#include "egs_rndm.h" #include #include @@ -149,6 +151,21 @@ class EGS_EXPORT EGS_BaseGeometry { */ virtual int isWhere(const EGS_Vector &x) = 0; + /* getNextGeom is the equivalent of getNextParticle but for the simulation + * object. Its goal is to determine the next state of the geometry, either + * by synchronizing itself to the source time parameter, or by sampling its + * own time parameter and updating itself accordingly if the source has + * provided no time index. + * + * This function has a non-empty implementation in 2 cases. + * + * 1) it is reimplemented in any composite geometry, where it will call next + * geom on all of its components + * + * 2) it is reimplemented in the dynamic geometry class. This is where the + * code will find the current (non static) state of the geometry. */ + virtual void getNextGeom(EGS_RandomGenerator *rndm) {}; + /*! \brief Find the bin to which \a xp belongs, given \a np bin edges \a p This static method uses a binary search and is provided for the @@ -679,7 +696,6 @@ class EGS_EXPORT EGS_BaseGeometry { return ++nref; }; - /*! \brief Decrease the reference count to this geometry Composite geometries should use this method to decrease @@ -731,6 +747,18 @@ class EGS_EXPORT EGS_BaseGeometry { /*! \brief Set the labels from an input string */ int setLabels(const string &inp); + virtual void updatePosition(EGS_Float time) { }; + + /* This method is essentially used to determine whether the simulation + * geometry contains a dynamic geometry. Like getNextGeom(), the only + * non-empty implementations of this function are in composite geometries + * (where it simply calls containsDynamic on its components), and in the + * dynamic geometry, where it will update the boolean reference to true and + * call on its base geometry. This function was conceived to be used in the + * view/viewcontrol (to determine whether time index objects are visible or + * hidden) */ + virtual void containsDynamic(bool &hasdynamic) { }; + protected: /*! \brief Number of local regions in this geometry @@ -763,7 +791,6 @@ class EGS_EXPORT EGS_BaseGeometry { */ int med; - /*! \brief Does this geometry have relative mass density scvaling? */ @@ -1116,7 +1143,6 @@ struct EGS_GeometryIntersections { \until make_depend That's all. - Here is the complete source code of the EGS_Box class.
The header file: \include geometry/egs_box/egs_box.h @@ -1124,7 +1150,6 @@ struct EGS_GeometryIntersections { \include geometry/egs_box/egs_box.cpp The Makefile: \include geometry/egs_box/Makefile - */ /* \example geometry/example1/geometry_example1.cpp @@ -1336,5 +1361,4 @@ struct EGS_GeometryIntersections { and the complete implementation: */ - #endif diff --git a/HEN_HOUSE/egs++/egs_base_source.h b/HEN_HOUSE/egs++/egs_base_source.h index 392eed3d8..5086f4408 100644 --- a/HEN_HOUSE/egs++/egs_base_source.h +++ b/HEN_HOUSE/egs++/egs_base_source.h @@ -24,6 +24,7 @@ # Author: Iwan Kawrakow, 2005 # # Contributors: Reid Townson +# Alexandre Demelo # ############################################################################### */ @@ -44,6 +45,7 @@ #include "egs_object_factory.h" #include "egs_functions.h" + #include #include #include "egs_math.h" @@ -85,7 +87,7 @@ class EGS_EXPORT EGS_BaseSource : public EGS_Object { * */ EGS_BaseSource(const string &Name="", EGS_ObjectFactory *f = 0) : - EGS_Object(Name,f) {}; + EGS_Object(Name,f), time_index(-1) {}; /*! \brief Construct a source from the input pointed to by \a inp. * @@ -96,7 +98,7 @@ class EGS_EXPORT EGS_BaseSource : public EGS_Object { * plus additional information as needed by the source being created. */ EGS_BaseSource(EGS_Input *input, EGS_ObjectFactory *f = 0) : - EGS_Object(input,f) {}; + EGS_Object(input,f), time_index(-1) {}; virtual ~EGS_BaseSource() {}; /*! \brief Get a short description of this source. @@ -184,9 +186,11 @@ class EGS_EXPORT EGS_BaseSource : public EGS_Object { * particle. Currently only makes sense for IAEA_PhspSource and * EGS_BeamSource. */ - virtual EGS_Float getMu() { - return -1; - }; + //virtual EGS_Float getTimeIndex() { + //return -1; + //}; + + //virtual void setTimeIndex(EGS_Float temp_time) {}; /*! \brief Store the source state into the stream \a data_out. * @@ -303,6 +307,32 @@ class EGS_EXPORT EGS_BaseSource : public EGS_Object { */ static void addKnownTypeId(const char *name); + /* Centralize the time index parameter so that it can be saved to and + * accessed from a single point. In almost any instance where a time index + * parameter is created, it is saved in the source object using + * setTimeIndex. When other objects would like to access the time index + * (most relevant example being the dynamic geometry checking if the source + * has provided a time index before setting its own), they can call + * getTimeIndex. In many cases, this method is indirectly called via the + * getTimeIndex in the application class, which returns the results of + * getTimeIndex call on the simulation source. Note there are two cases + * which behave slightly differently and may need some modifications. While + * the dynamicSource time index implementation was completely absorbed into + * the basesource, this was not done for the beam source and the iaea_phsp + * source, as they had slightly more involved time index implementations + * that seemed best left alone. They do not have setTimeIndex methods, \ + * and may not return the right thing if we set the time using the geometry, + * as calling get may try to get the beam or iaea_phsp time index and not + * the basesource index we set with the geometry. */ + EGS_Float getTimeIndex() { + return time_index; + + }; + + void setTimeIndex(EGS_Float temp_time) { + time_index=temp_time; + }; + protected: /*! \brief A short source description. @@ -312,6 +342,11 @@ class EGS_EXPORT EGS_BaseSource : public EGS_Object { */ string description; + /*! \brief time index corresponding to a particle. This stores the current + * time index for all objects in the simulation (with the potential + * exception of beam and iaea_phsp source) */ + EGS_Float time_index; + }; /*! \brief Base class for energy spectra. All energy spectra in the EGSnrc @@ -567,7 +602,7 @@ class EGS_EXPORT EGS_BaseSimpleSource : public EGS_BaseSource { */ EGS_BaseSimpleSource(int Q, EGS_BaseSpectrum *Spec, const string &Name="", EGS_ObjectFactory *f=0) : - EGS_BaseSource(Name,f), q(Q), s(Spec), count(0) { }; + EGS_BaseSource(Name,f), q(Q), s(Spec), count(0) {}; /*! \brief Construct a 'simple' particle source from the information * pointed to by \a input. @@ -777,6 +812,7 @@ class EGS_EXPORT EGS_BaseSimpleSource : public EGS_BaseSource { return true; }; + protected: /*! Set the particle latch. @@ -801,6 +837,8 @@ class EGS_EXPORT EGS_BaseSimpleSource : public EGS_BaseSource { /*! \brief Number of statistically independent particles delivered so far.*/ EGS_I64 count; + + }; /*! \brief A template source creation function. diff --git a/HEN_HOUSE/egs++/egs_particle_track.cpp b/HEN_HOUSE/egs++/egs_particle_track.cpp index 7ad165aa2..d8a0db489 100644 --- a/HEN_HOUSE/egs++/egs_particle_track.cpp +++ b/HEN_HOUSE/egs++/egs_particle_track.cpp @@ -24,6 +24,7 @@ # Author: Georgi Gerganov, 2009 # # Contributors: Iwan Kawrakow +# Alexandre Demelo # ############################################################################### */ @@ -35,6 +36,10 @@ #include "egs_particle_track.h" +#include +#include + + void EGS_ParticleTrack::grow() { // calculate the new size of the vertex array @@ -70,13 +75,20 @@ void EGS_ParticleTrack::clearTrack() { m_nVertices = 0; } -int EGS_ParticleTrack::writeTrack(ofstream *trsp) { +int EGS_ParticleTrack::writeTrack(ofstream *trsp, bool incltime) { // no need to write the track if it has less than 2 vertices ... if (m_nVertices < 2) { return 1; } + trsp->write((char *)&m_nVertices, sizeof(int)); trsp->write((char *)m_pInfo, sizeof(ParticleInfo)); + + // If user indicates to include time in inputfile then write the time index + // to the tracks file after the particle info. + if (incltime) { + trsp->write((char *) &time_index, sizeof(EGS_Float)); + } for (int i = 0; i < m_nVertices; i++) { trsp->write((char *)m_track[i],sizeof(Vertex)); } @@ -99,7 +111,6 @@ void EGS_ParticleTrack::addVertex(Vertex *x) { } } -///////////////////////// EGS_ParticleTrackContainer //////////////////////// void EGS_ParticleTrackContainer::startNewTrack() { // buffer full ? flush it @@ -139,7 +150,7 @@ void EGS_ParticleTrackContainer::flushBuffer() { // if the particle is not being scored anymore if (!m_isScoring[i]) { // output it to the file and free memory - if (!m_buffer[i]->writeTrack(m_trspFile)) { + if (!m_buffer[i]->writeTrack(m_trspFile,incltime)) { m_totalTracks++; } m_buffer[i]->clearTrack(); @@ -204,6 +215,7 @@ int EGS_ParticleTrackContainer::readDataFile(const char *filename) { func_name, filename); return -1; } + data->read((char *)&m_totalTracks, sizeof(int)); egsInformation("%s: Reading %d tracks from '%s' ...\n", func_name, m_totalTracks, filename); m_nTracks = 0; @@ -214,6 +226,7 @@ int EGS_ParticleTrackContainer::readDataFile(const char *filename) { m_isScoring = new bool [m_totalTracks]; for (int i = 0; i < m_totalTracks; i++) { int nvertices; + EGS_Float time; m_buffer[i] = new EGS_ParticleTrack(); m_stackMap[i] = -1; m_isScoring[i] = false; @@ -222,6 +235,12 @@ int EGS_ParticleTrackContainer::readDataFile(const char *filename) { totalVertices += nvertices; data->read((char *)pinfo,sizeof(EGS_ParticleTrack::ParticleInfo)); startNewTrack(pinfo); + + // If user indicates to include time in inputfile then read the time + // index from the tracks file after the particle info. + if (incltime) { + data->read((char *)&time,sizeof(EGS_Float)); + } for (int j = 0; j < nvertices; j++) { EGS_ParticleTrack::Vertex *v = new EGS_ParticleTrack::Vertex(); data->read((char *)v,sizeof(EGS_ParticleTrack::Vertex)); @@ -266,4 +285,87 @@ void EGS_ParticleTrackContainer::reportResults(bool with_header) { } } egsInformation(" Output file name: %s\n\n", m_trspFilename.c_str()); + if (incltime) { + tracksFileSort(); + } + readDataFile(m_trspFilename.c_str()); +} + +void EGS_ParticleTrackContainer::tracksFileSort() { + + // Function to sort tracks file entries by time index. This is done by + // reading time indices, saving their positions in a vector of pairs, + // sorting the vector by time index, and then rewriting tracks into a new + // sorted file. + + const char *func_name = "EGS_ParticleTrackContainer::tracksFileSort()"; + const char *trackfile = m_trspFilename.c_str(); + ifstream *data = new ifstream(trackfile, ios::binary); + + // New sorted tracks file where data from the original tracksfile will be + // rewritten. + string outstring = "sorted_trackfile.syncptracks"; + const char *outname = outstring.c_str(); + ofstream *sortout = new ofstream(outname, ios::binary); + + int totalTrackNum; + data->read((char *)&totalTrackNum, sizeof(int)); + egsInformation("%s: Sorting %d tracks from '%s' by time index ...\n", func_name, totalTrackNum, trackfile); + sortout->write((char *)&totalTrackNum, sizeof(int)); + + // Defining vector of pairs which will be used to sort. + vector> vect; + data->seekg(sizeof(int), ios_base::beg); // Skip the total number of tracks in the beginning. + + for (int i = 0; i < m_totalTracks; i++) { + int nvertices; + EGS_ParticleTrack::ParticleInfo *pinfo = new EGS_ParticleTrack::ParticleInfo(0); + EGS_Float time; + int position = data->tellg(); // Get the position in the file (right before a new track starts). + + // Read in usual track info. + data->read((char *)&nvertices, sizeof(int)); + data->read((char *)pinfo, sizeof(EGS_ParticleTrack::ParticleInfo)); + data->read((char *)&time, sizeof(EGS_Float)); + vect.push_back(make_pair(time, position)); // Add the time index/file position pair to the vector. + + // Determine the size of all the vertices for track i and skip these to + // get to the next track. + int vertsize = nvertices * sizeof(EGS_ParticleTrack::Vertex); + data->seekg(vertsize, ios_base::cur); + } + + // Using the built-in C++ sort function, sort all the pairs by their time + // index (first element in the pair). + sort(vect.begin(), vect.end()); + + // Now, write the sorted tracks file by jumping around the unsorted file + // using the positions in the pairs vector. + for (int i = 0; i < m_totalTracks; i++) { + int nvertices; + EGS_ParticleTrack::ParticleInfo *pinfo = new EGS_ParticleTrack::ParticleInfo(0); + EGS_Float time; + data->seekg(vect[i].second); // Go to file position right before track i. + + // Read in all track info for track i from unsorted. + data->read((char *)&nvertices, sizeof(int)); + data->read((char *)pinfo, sizeof(EGS_ParticleTrack::ParticleInfo)); + data->read((char *)&time, sizeof(EGS_Float)); + + // Write all read track info to the sorted file (at the end of the file). + sortout->write((char *)&nvertices, sizeof(int)); + sortout->write((char *)pinfo, sizeof(EGS_ParticleTrack::ParticleInfo)); + sortout->write((char *)&time, sizeof(EGS_Float)); + + for (int j = 0; j < nvertices; j++) { + // Read the nvertices from the unsorted file and write them to the end of the sorted file. + EGS_ParticleTrack::Vertex *v = new EGS_ParticleTrack::Vertex(); + data->read((char *)v, sizeof(EGS_ParticleTrack::Vertex)); + sortout->write((char *)v, sizeof(EGS_ParticleTrack::Vertex)); + } + } + + sortout->close(); + int removal = remove(trackfile); // Delete unsorted file. + int renaming = rename(outname, trackfile); // Rename sorted file to the unsorted file's old name. } diff --git a/HEN_HOUSE/egs++/egs_particle_track.h b/HEN_HOUSE/egs++/egs_particle_track.h index acfb2f100..47c23022d 100644 --- a/HEN_HOUSE/egs++/egs_particle_track.h +++ b/HEN_HOUSE/egs++/egs_particle_track.h @@ -24,6 +24,7 @@ # Author: Georgi Gerganov, 2009 # # Contributors: Iwan Kawrakow +# Alexandre Demelo # ############################################################################### */ @@ -116,7 +117,7 @@ class EGS_EXPORT EGS_ParticleTrack { Writes the particle information block m_pInfo, the number of vertices in this track and finally all vertex data. */ - int writeTrack(ofstream *trsp); + int writeTrack(ofstream *trsp, bool incltime); /*! \brief Add additional point of interaction (Vertex) to the track. */ void addVertex(Vertex *x); @@ -143,12 +144,22 @@ class EGS_EXPORT EGS_ParticleTrack { return m_pInfo; }; + void setTimeIndex(EGS_Float time) { + time_index = time; + }; + + /*! \brief Get the time index of the particle being tracked. */ + EGS_Float getTimeIndex() { + return time_index; + }; + protected: /*! \brief Resize the array containing the vertices. */ void grow(); ParticleInfo *m_pInfo; //!< type of the tracked particle + EGS_Float time_index; int m_size; //!< current size of the vertex array int m_nVertices; //!< current number of vertices in track @@ -181,12 +192,15 @@ class EGS_EXPORT EGS_ParticleTrackContainer { \a buf_size which defines how many tracks the container will store before flushing them to the output file. */ - EGS_ParticleTrackContainer(const char *fname, int buf_size) : m_nEvents(0), + EGS_ParticleTrackContainer(const char *fname, int buf_size, bool time_bool) : m_nEvents(0), m_nTracks(0), m_totalTracks(0), m_isScoring(NULL), m_bufferSize(0), m_buffer(0), m_trspFile(NULL) { m_bufferSize = buf_size; - // initialize the arrays + // Initialize the arrays. The incltime variable set from the time_bool + // passed from the constructor (track scoring object gets incltime from + // input file and calls and passes) + incltime = time_bool; m_buffer = new EGS_ParticleTrack* [m_bufferSize]; m_stackMap = new int [m_bufferSize]; m_isScoring = new bool [m_bufferSize]; @@ -202,6 +216,7 @@ class EGS_EXPORT EGS_ParticleTrackContainer { int dummy = 0; // at the end this will be replaced with the number of recorded tracks m_trspFile->write((char *)&dummy, sizeof(int)); + }; /*! \brief The Destructor. Deallocate all allocated memory. */ @@ -317,6 +332,10 @@ class EGS_EXPORT EGS_ParticleTrackContainer { m_buffer[m_nTracks-1]->setParticleInfo(p); } + void setCurrentTimeIndex(EGS_Float time) { + m_buffer[m_nTracks-1]->setTimeIndex(time); + } + /*! \brief Load particle data from the file called \a filename .*/ int readDataFile(const char *filename); @@ -325,6 +344,8 @@ class EGS_EXPORT EGS_ParticleTrackContainer { */ void reportResults(bool with_header = true); + void tracksFileSort(); + protected: /*! \brief Write all track data to the file. @@ -341,6 +362,7 @@ class EGS_EXPORT EGS_ParticleTrackContainer { */ void updateHeader(); + bool incltime; int m_nEvents; //!< number of events scored int m_nTracks; //!< number of tracks currently in memory int m_totalTracks; //!< total number of tracks registered diff --git a/HEN_HOUSE/egs++/egs_shapes.cpp b/HEN_HOUSE/egs++/egs_shapes.cpp index d69356e0a..1a6b034ce 100644 --- a/HEN_HOUSE/egs++/egs_shapes.cpp +++ b/HEN_HOUSE/egs++/egs_shapes.cpp @@ -56,6 +56,7 @@ EGS_BaseShape *EGS_BaseShape::createShape(EGS_Input *i) { shape_creator.addKnownObject(new EGS_CylinderShape); } __shape_count++; + EGS_Object *o = shape_creator.createSingleObject(i,"createShape",true); return dynamic_cast(o); } @@ -189,11 +190,11 @@ EGS_Object *EGS_CylinderShape::createObject(EGS_Input *input) { bool has_A = (err2 == 0 && A.size() == 3); EGS_CylinderShape *result; if (has_Xo && has_A) result = new EGS_CylinderShape(r,H, - EGS_Vector(Xo[0],Xo[1],Xo[2]),EGS_Vector(A[0],A[1],A[2])); + EGS_Vector(Xo[0],Xo[1],Xo[2]),EGS_Vector(A[0],A[1],A[2])); else if (has_Xo) result = new EGS_CylinderShape(r,H, - EGS_Vector(Xo[0],Xo[1],Xo[2])); + EGS_Vector(Xo[0],Xo[1],Xo[2])); else if (has_A) result = new EGS_CylinderShape(r,H, - EGS_Vector(0,0,0),EGS_Vector(A[0],A[1],A[2])); + EGS_Vector(0,0,0),EGS_Vector(A[0],A[1],A[2])); else { result = new EGS_CylinderShape(r,H); } diff --git a/HEN_HOUSE/egs++/egs_shapes.h b/HEN_HOUSE/egs++/egs_shapes.h index dbefb21b0..a5b164046 100644 --- a/HEN_HOUSE/egs++/egs_shapes.h +++ b/HEN_HOUSE/egs++/egs_shapes.h @@ -205,6 +205,14 @@ class EGS_EXPORT EGS_BaseShape : public EGS_Object { return false; }; + /* getNextShapePosition is the equivalent of getNextParticle but for a shape object. Its goal is to determine the next state of the geometry, either by synchronizing itself to the + * source time parameter, or by sampling it's own time parameter and updating itself accordingly if the source has provided no time index. + * + *This function has a non-empty implementation in 2 cases. + *1) it is re implemented in any composite shape, where it will call getNextShapePosition on all of its components + *2) it is re implemented in the dynamic shape class. This is where the code will find the current (non static) state of the shape. */ + virtual void getNextShapePosition(EGS_RandomGenerator *rndm) {}; + /*! Get a random direction given a source position \a xo. * * This method is used to sample random directions by picking a @@ -233,6 +241,10 @@ class EGS_EXPORT EGS_BaseShape : public EGS_Object { return 1; }; + /*! \brief Update the position of the shape if it is in motion + */ + virtual void updatePosition(EGS_Float time) { }; + protected: EGS_AffineTransform *T; //!< The affine transformation attached to the shape @@ -388,7 +400,7 @@ class EGS_EXPORT EGS_BoxShape : public EGS_BaseShape { /*! \brief Destructor. Does nothing */ ~EGS_BoxShape() { }; - /*! \brief Returns a point uniformely distributed within the box. */ + /*! \brief Returns a point uniformly distributed within the box. */ EGS_Vector getPoint(EGS_RandomGenerator *rndm) { EGS_Vector v(ax*(rndm->getUniform()-0.5), ay*(rndm->getUniform()-0.5), @@ -408,7 +420,7 @@ class EGS_EXPORT EGS_BoxShape : public EGS_BaseShape { return true; }; - /*! \brief Sets the direction \a u by picking a random point uniformely the + /*! \brief Sets the direction \a u by picking a random point uniformly the * on the box surface. * \sa EGS_BaseShape::getPointSourceDirection() */ @@ -538,7 +550,7 @@ class EGS_EXPORT EGS_SphereShape : public EGS_BaseShape { return true; }; - /*! \brief Sets the direction \a u by picking a random point uniformely + /*! \brief Sets the direction \a u by picking a random point uniformly * on the sphere surface. * \sa EGS_BaseShape::getPointSourceDirection() */ diff --git a/HEN_HOUSE/egs++/egs_transformations.cpp b/HEN_HOUSE/egs++/egs_transformations.cpp index 23b8c0d9c..0eb590d73 100644 --- a/HEN_HOUSE/egs++/egs_transformations.cpp +++ b/HEN_HOUSE/egs++/egs_transformations.cpp @@ -23,7 +23,7 @@ # # Author: Iwan Kawrakow, 2005 # -# Contributors: +# Contributors: Alexandre Demelo # ############################################################################### */ @@ -105,3 +105,39 @@ EGS_AffineTransform *EGS_AffineTransform::getTransformation(EGS_Input *i) { } return result; } + +// Overload for dynamic geometry +EGS_AffineTransform *EGS_AffineTransform::getTransformation(vector translation, vector rotation) { + + /* The following method is used by the dynamic geometry class to create a + * transformation corresponding to a certain translation and rotation + * vector. These vectors are determined through the control points and the + * sampled time value. The sampled transformation coordinates are passed + * here to build a transformation */ + + // First check that appropriate values have been provided. May work with + // either 2 or 3 rotation parameters, and exactly 3 translation parameters + if (translation.size()!=3 || (rotation.size()!=2 && rotation.size()!=3)) { + egsWarning("getTransformation: invalid transformation parameters\n"); + return 0; + } + + EGS_AffineTransform *result; //the returned transformation + + // EGS_vector holding translation parameters is defined to pass to the + // EGS_AffineTransform constructor + EGS_Vector t = EGS_Vector(translation[0],translation[1],translation[2]); + + // Check the size of the rotation vector provided and call the appropriate + // EGS_AffineTransform constructor (converting rotation vector to matrix) + if (rotation.size() == 2) { + result = new EGS_AffineTransform(EGS_RotationMatrix(rotation[0],rotation[1]),t); + } + else if (rotation.size() == 3) { + result = new EGS_AffineTransform(EGS_RotationMatrix(rotation[0],rotation[1],rotation[2]),t); + } + else { + result = new EGS_AffineTransform(EGS_RotationMatrix(),t); + } + return result; +} diff --git a/HEN_HOUSE/egs++/egs_transformations.h b/HEN_HOUSE/egs++/egs_transformations.h index dac540615..5c91321d1 100644 --- a/HEN_HOUSE/egs++/egs_transformations.h +++ b/HEN_HOUSE/egs++/egs_transformations.h @@ -24,6 +24,7 @@ # Author: Iwan Kawrakow, 2005 # # Contributors: +# Alexandre Demelo # ############################################################################### */ @@ -44,7 +45,9 @@ #include "egs_math.h" #include "egs_functions.h" + #include +#include using namespace std; @@ -652,6 +655,7 @@ class EGS_EXPORT EGS_AffineTransform { */ static EGS_AffineTransform *getTransformation(EGS_Input *inp); + static EGS_AffineTransform *getTransformation(vector trnsl, vector rot); }; #endif diff --git a/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.cpp b/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.cpp index c46d9d501..cb9d8dc38 100644 --- a/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.cpp +++ b/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.cpp @@ -27,6 +27,7 @@ # Contributors: Marc Chamberland # Rowan Thomson # Dave Rogers +# Alexandre Demelo # ############################################################################### # @@ -606,6 +607,38 @@ int EGS_AEnvelope::getMaxStep() const { return nstep + inscribed_geoms.size(); }; +void EGS_AEnvelope::getNextGeom(EGS_RandomGenerator *rndm) { + // calls getNextGeom on its component geometries to update dynamic + // geometries in the simulation + for (int j=0; j< inscribed_geoms.size(); ++j) { + inscribed_geoms[j]->getNextGeom(rndm); + } + base_geom->getNextGeom(rndm); +}; + +void EGS_AEnvelope::updatePosition(EGS_Float time) { + // calls updatePosition on its component geometries to update dynamic + // geometries in the simulation + for (int j=0; j< inscribed_geoms.size(); j++) { + inscribed_geoms[j]->updatePosition(time); + } + base_geom->updatePosition(time); +}; + +void EGS_AEnvelope::containsDynamic(bool &hasdynamic) { + // calls containsDynamic on its component geometries (only calls if + // hasDynamic is false, as if it is true we already found one) + for (int j=0; j< inscribed_geoms.size(); j++) { + if (!hasdynamic) { + inscribed_geoms[j]->containsDynamic(hasdynamic); + } + } + if (!hasdynamic) { + base_geom->containsDynamic(hasdynamic); + } +}; + + EGS_Float EGS_AEnvelope::getVolume(int ireg) { if (ireg < 0) { @@ -892,9 +925,9 @@ void EGS_ASwitchedEnvelope::activateByIndex(int inscribed_index) { void EGS_ASwitchedEnvelope::deactivateByIndex(int inscribed_index) { vector::iterator loc = find( - active_inscribed.begin(), active_inscribed.end(), - inscribed_geoms[inscribed_index] - ); + active_inscribed.begin(), active_inscribed.end(), + inscribed_geoms[inscribed_index] + ); if (loc != active_inscribed.end()) { active_inscribed.erase(loc); } diff --git a/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.h b/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.h index 4a826adc5..bee4ed12d 100644 --- a/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.h +++ b/HEN_HOUSE/egs++/geometry/egs_autoenvelope/egs_autoenvelope.h @@ -27,6 +27,7 @@ # Contributors: Marc Chamberland # Rowan Thomson # Dave Rogers +# Alexandre Demelo # ############################################################################### # @@ -402,6 +403,12 @@ class EGS_AENVELOPE_EXPORT EGS_AEnvelope : public EGS_BaseGeometry { int getMaxStep() const; + void getNextGeom(EGS_RandomGenerator *rndm); + + void updatePosition(EGS_Float time); + + void containsDynamic(bool &hasdynamic); + virtual EGS_Float getVolume(int ireg); virtual EGS_Float getCorrectionRatio(int ireg); diff --git a/HEN_HOUSE/egs++/geometry/egs_cd_geometry/egs_cd_geometry.h b/HEN_HOUSE/egs++/geometry/egs_cd_geometry/egs_cd_geometry.h index 228ef2975..6cd50d6cc 100644 --- a/HEN_HOUSE/egs++/geometry/egs_cd_geometry/egs_cd_geometry.h +++ b/HEN_HOUSE/egs++/geometry/egs_cd_geometry/egs_cd_geometry.h @@ -26,6 +26,7 @@ # Contributors: Ernesto Mainegra-Hing # Frederic Tessier # Reid Townson +# Alexandre Demelo # ############################################################################### */ @@ -819,6 +820,32 @@ class EGS_CDGEOMETRY_EXPORT EGS_CDGeometry : public EGS_BaseGeometry { return nstep + 1; } + void getNextGeom(EGS_RandomGenerator *rndm) { //calls getNextGeom on its component geometries to update dynamic geometries in the simulation + for (int j=0; jregions(); ++j) { + g[j]->getNextGeom(rndm); + } + bg->getNextGeom(rndm); + }; + + void updatePosition(EGS_Float time) {//calls updatePosition on its component geometries to update dynamic geometries in the simulation + for (int j=0; jregions(); j++) { + g[j]->updatePosition(time); + } + bg->updatePosition(time); + }; + + void containsDynamic(bool &hasdynamic) {//calls containsDynamic on its component geometries (only calls if hasDynamic is false, as if it is true we already found one) + for (int j=0; jregions(); j++) { + if (!hasdynamic) { + g[j]->containsDynamic(hasdynamic); + } + } + if (!hasdynamic) { + bg->containsDynamic(hasdynamic); + } + }; + + bool hasBooleanProperty(int ireg, EGS_BPType prop) const { if (ireg >= 0 && ireg < nreg) { int ibase = ireg/nmax; diff --git a/HEN_HOUSE/egs++/geometry/egs_conez/egs_conez.h b/HEN_HOUSE/egs++/geometry/egs_conez/egs_conez.h index d72d9ccc6..a85016a68 100644 --- a/HEN_HOUSE/egs++/geometry/egs_conez/egs_conez.h +++ b/HEN_HOUSE/egs++/geometry/egs_conez/egs_conez.h @@ -63,7 +63,7 @@ template class EGS_CONEZ_EXPORT EGS_ConezT : public EGS_BaseGeometry { protected: - EGS_Float *theta,*cos_t,*cos2_t,*sin_t,*sin2_t; + EGS_Float *theta, *cos_t, *cos2_t, *sin_t, *sin2_t; // for conc-cones, all apices coincide EGS_Vector xo; // projection operator diff --git a/HEN_HOUSE/egs++/geometry/egs_dynamic_geometry/Makefile b/HEN_HOUSE/egs++/geometry/egs_dynamic_geometry/Makefile new file mode 100644 index 000000000..d12b11b06 --- /dev/null +++ b/HEN_HOUSE/egs++/geometry/egs_dynamic_geometry/Makefile @@ -0,0 +1,44 @@ + +############################################################################### +# +# EGSnrc egs++ makefile to build dynamic_geometry geometry +# 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: Alexandre Demelo, 2023 +# +# Contributors: +# +############################################################################### + + +include $(EGS_CONFIG) +include $(SPEC_DIR)egspp.spec +include $(SPEC_DIR)egspp_$(my_machine).conf + +DEFS = $(DEF1) -DBUILD_DYNAMIC_GEOMETRY_DLL + +library = egs_dynamic_geometry +lib_files = egs_dynamic_geometry +my_deps = egs_transformations.h +extra_dep = $(addprefix $(DSOLIBS), $(my_deps)) + +include $(SPEC_DIR)egspp_libs.spec + +$(make_depend) diff --git a/HEN_HOUSE/egs++/geometry/egs_dynamic_geometry/egs_dynamic_geometry.cpp b/HEN_HOUSE/egs++/geometry/egs_dynamic_geometry/egs_dynamic_geometry.cpp new file mode 100644 index 000000000..822cd98a8 --- /dev/null +++ b/HEN_HOUSE/egs++/geometry/egs_dynamic_geometry/egs_dynamic_geometry.cpp @@ -0,0 +1,629 @@ +/* +############################################################################### +# +# EGSnrc egs++ dynamic_geometry geometry +# 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: Alexandre Demelo, 2023 +# +# Contributors: Reid Townson +# +############################################################################### +*/ + + +/*! \file egs_dynamic_geometry.cpp + * \brief A dynamic geometry + * \RT + */ + +#include "egs_dynamic_geometry.h" +#include "egs_input.h" +#include "egs_functions.h" + +// ---------------------------------------------------------------------------- +// Implementation of EGS_DynamicGeometry methods + +// ---------------------------------------------------------------------------- +// Avoid using deprecated methods from the base class +void EGS_DynamicGeometry::setMedia(EGS_Input *, int, const int *) { + egsWarning("EGS_DynamicGeometry::setMedia: don't use this method. Use the\n" + " setMedia() methods of the geometry objects that make up this geometry\n"); +} + +void EGS_DynamicGeometry::setRelativeRho(int start, int end, EGS_Float rho) { + setRelativeRho(0); +} + +void EGS_DynamicGeometry::setRelativeRho(EGS_Input *) { + egsWarning("EGS_DynamicGeometry::setRelativeRho(): don't use this " + "method. Use the \n setRelativeRho() methods of the underlying " + "geometry\n"); +} + +void EGS_DynamicGeometry::setBScaling(int start, int end, EGS_Float rho) { + setBScaling(0); +} + +void EGS_DynamicGeometry::setBScaling(EGS_Input *) { + egsWarning("EGS_DynamicGeometry::setBScaling(): don't use this " + "method. Use the \n setBScaling() methods of the underlying " + "geometry\n"); +} + +// ---------------------------------------------------------------------------- +// External C function to create EGS_DynamicGeometry +extern "C" { + + /*! \brief Create a dynamic geometry from input. + * + * This function is part of the EGSnrc dynamic geometry plugin interface. + * \param input The input specifying the dynamic geometry. + * \return A pointer to the created EGS_DynamicGeometry object. + */ + EGS_DYNAMIC_GEOMETRY_EXPORT EGS_BaseGeometry *createGeometry(EGS_Input *input) { + EGS_BaseGeometry *g = 0; + EGS_Input *ij = input->takeInputItem("geometry", false); + if (ij) { + g = EGS_BaseGeometry::createSingleGeometry(ij); + delete ij; + if (!g) { + egsWarning("createGeometry(EGS_DynamicGeometry): got a null pointer as a geometry?\n"); + return 0; + } + } + if (!g) { + string gname; + int err = input->getInput("my geometry", gname); + if (err) { + egsWarning("createGeometry(EGS_DynamicGeometry): my geometry must be defined\n either inline or using 'my geometry = some_name'\n"); + return 0; + } + g = EGS_BaseGeometry::getGeometry(gname); + if (!g) { + egsWarning("createGeometry(EGS_DynamicGeometry): no geometry named %s is defined\n", gname.c_str()); + return 0; + } + } + g->ref(); + + // Now getting motion information for dynamic component + EGS_Input *dyninp = input->takeInputItem("motion"); + + if (dyninp) { + EGS_DynamicGeometry *result = new EGS_DynamicGeometry(g, dyninp); + result->setName(input); + result->setBoundaryTolerance(input); + result->setLabels(input); + return result; + } + else { + egsWarning("EGS_DynamicGeometry: no control points input.\n"); + return 0; + } + } + + /*! \brief Get label regions for dynamic geometry. + * + * This function is part of the EGSnrc dynamic geometry plugin interface. + * \param str Label string. + * \param regs Vector to store label regions. + */ + void EGS_DynamicGeometry::getLabelRegions(const string &str, vector ®s) { + // label defined in the geometry being transformed + g->getLabelRegions(str, regs); + + // label defined in self (transformation input block) + EGS_BaseGeometry::getLabelRegions(str, regs); + } + + /*! \brief Build dynamic geometry from input. + * + * This function is part of the EGSnrc dynamic geometry plugin interface. + * \param g Pointer to the base geometry. + * \param dyninp Input specifying the dynamic geometry. + */ + void EGS_DynamicGeometry::buildDynamicGeometry(EGS_BaseGeometry *g, EGS_Input *dyninp) { + stringstream itos; + ncpts = 0; + vector point; + EGS_ControlPoint cpt; + int err; + int icpts = 1; + itos << icpts; + string inputTag = "control point"; + string inputTag_backCompat = "control point " + itos.str(); + EGS_Input *currentInput; + int rotsize = 0; + + // Control points read one by one from motion block in dynamic geometry + // definition, and saved as a vector to points + while (true) { + currentInput = dyninp->takeInputItem(inputTag); + if (!currentInput || currentInput->getInput(inputTag, point)) { + currentInput = dyninp->takeInputItem(inputTag_backCompat); + if (!currentInput || currentInput->getInput(inputTag_backCompat, point)) { + delete currentInput; + break; + } + } + delete currentInput; + + // Checking the size to make sure it is a valid control point input + if (point.size() != 6 && point.size() != 7) { + egsWarning("EGS_DynamicGeometry: Control point %i must specify either 6 or 7 values.\n", icpts); + } + else { + if (ncpts == 0) { + // Variable to make sure all control point definitions have consistent formats + rotsize = point.size(); + } + + // Make sure each time index is larger than the last + if (ncpts > 0 && point[0] < cpts[ncpts - 1].time) { + egsWarning("EGS_DynamicGeometry: Time index of control point %i < time index of control point %i\n", icpts, ncpts); + } + + // Checks that time index is valid (larger than zero) + else if (point[0] < 0.0) { + egsWarning("EGS_DynamicGeometry: Time index of control point %i < 0.0\n", icpts); + } + + // Checks that control point formats follow the first + else if (ncpts > 0 && point.size() != rotsize) { + egsWarning("EGS_DynamicGeometry: Rotation definition inconsistent \n"); + } + else { + ncpts++; + if (ncpts == 1 && point[0] > 0.0) { + egsWarning("EGS_DynamicGeometry: Time index of control point 1 > 0.0. This will generate many warning messages.\n"); + } + vector T_vect; // Vector storing translation information + vector R_vect; // Vector storing rotation information + + // Add translation coordinates to the translation vector + T_vect.push_back(point[1]); + T_vect.push_back(point[2]); + T_vect.push_back(point[3]); + + // Add rotation coordinates to rotation vector (cases + // differentiate cpt formats, 6 is 2 rotation parameters, 7 + // is 3) In each case vector order is determined by format + // of EGS_rotationMatrix constructor (EGS_Float arguments) + if (point.size() == 6) { + R_vect.push_back(point[6]); // Rotation about z + R_vect.push_back(point[4]); // Rotation about x + } // 2 number case + if (point.size() == 7) { + R_vect.push_back(point[4]); + R_vect.push_back(point[5]); // Rotation about y + R_vect.push_back(point[6]); // Rotation about z + } // 3 number case + + // Add it to the vector of control points + cpt.time = point[0]; + cpt.trnsl = T_vect; + cpt.rot = R_vect; + cpts.push_back(cpt); + icpts++; + itos.str(""); + itos << icpts; + + // Define next control point i string for getInput in while condition + inputTag_backCompat = "control point " + itos.str(); + } + } + } + if (ncpts <= 1) { + egsFatal("EGS_DynamicGeometry: Not enough or missing control points.\n"); + } + if (cpts[ncpts - 1].time == 0.0) { + egsFatal("EGS_DynamicGeometry: Time index of last control point = 0. Something's wrong.\n"); + } + else { + + // Normalize time index to max. value + for (int i = 0; i <= ncpts - 1; i++) { + cpts[i].time /= cpts[ncpts - 1].time; + } + } + + // Sets position to initial time in egs_view upon opening + updatePosition(0); + } + + /*! \brief Get interpolated coordinates based on the time index + * + * This function is part of the EGSnrc dynamic geometry plugin interface. + * \param rand Random number for sampling. + * \param gipt Control point struct used to track the current state of the geometry. + * \return 0 if successful, 1 otherwise. + */ + int EGS_DynamicGeometry::getCoordGeom(EGS_Float rand, EGS_ControlPoint &gipt) { + int iindex = 0; + int i; + + // The following loop determines which 2 control points the current time + // index falls between + for (i = 0; i < ncpts; i++) { + if (rand < cpts[i].time - epsilon) { + iindex = i; + break; + } + } + + if (i == ncpts) { + egsWarning("EGS_DynamicGeometry: could not locate control point.\n"); + return 1; + } + + // Below 3 vectors are defined, a vector containing the lower bound + // translation coordinates, a vector containing the upper bound + // translation coordinates, and a vector for the sampled translation + // coordinates + vector translation_LB = cpts[iindex - 1].trnsl; + vector translation_UB = cpts[iindex].trnsl; + vector translation_samp; + + // The following is a factor (between 0 and 1) used for sampling. + // Essentially, it represents the fractional position within the + // interval + EGS_Float factor = (rand - cpts[iindex - 1].time) / (cpts[iindex].time - cpts[iindex - 1].time); + + // Translations are given as the lower bound plus the length of the + // interval multiplied by the fractional position factor. So its + // essentially lower bound + n% of the interval length + translation_samp.push_back(translation_LB[0] + (translation_UB[0] - translation_LB[0]) * factor); + translation_samp.push_back(translation_LB[1] + (translation_UB[1] - translation_LB[1]) * factor); + translation_samp.push_back(translation_LB[2] + (translation_UB[2] - translation_LB[2]) * factor); + + // Update the translation coordinates in the current state control point object + gipt.trnsl = translation_samp; + + // Again 3 vectors defined, lower bound, upper bound, and sampled, this + // time for rotation coordinates + vector rotation_LB = cpts[iindex - 1].rot; + vector rotation_UB = cpts[iindex].rot; + vector rotation_samp; + + // Now set rotations. Current coordinates computed as before + // (lowerbound+n% of interval length), but now we also convert degrees + // to radians These must be done case by case (since order of arguments + // matters) + if (cpts[iindex].rot.size() == 2) { + rotation_samp.push_back((rotation_LB[0] + (rotation_UB[0] - rotation_LB[0]) * factor) * (M_PI / 180)); + rotation_samp.push_back((rotation_LB[1] + (rotation_UB[1] - rotation_LB[1]) * factor) * (M_PI / 180)); + } + else if (cpts[iindex].rot.size() == 3) { + rotation_samp.push_back((rotation_LB[0] + (rotation_UB[0] - rotation_LB[0]) * factor) * (M_PI / 180)); + rotation_samp.push_back((rotation_LB[1] + (rotation_UB[1] - rotation_LB[1]) * factor) * (M_PI / 180)); + rotation_samp.push_back((rotation_LB[2] + (rotation_UB[2] - rotation_LB[2]) * factor) * (M_PI / 180)); + } + else { + egsWarning("EGS_DynamicGeometry: Invalid number of rotation parameters\n"); + } + + // Update the rotation coordinates in the current state control point object + gipt.rot = rotation_samp; + + return 0; + } + + /*! \brief Compute intersections for dynamic geometry. + * + * This function is part of the EGSnrc dynamic geometry plugin interface. + * \param ireg Current region index. + * \param n Number of particles. + * \param x Position vector. + * \param u Direction vector. + * \param isections Pointer to EGS_GeometryIntersections object. + * \return Number of intersections. + */ + int EGS_DynamicGeometry::computeIntersections(int ireg, int n, const EGS_Vector &x, + const EGS_Vector &u, EGS_GeometryIntersections *isections) { + EGS_Vector xt(x), ut(u); + T.inverseTransform(xt); + T.rotateInverse(ut); + return g->computeIntersections(ireg, n, xt, ut, isections); + } + + /** + * @brief Check if the given region index corresponds to a real region in the dynamic geometry. + * + * @param ireg The region index. + * @return True if the region is real, false otherwise. + */ + bool EGS_DynamicGeometry::isRealRegion(int ireg) const { + return g->isRealRegion(ireg); + } + + /** + * @brief Check if a point is inside the dynamic geometry. + * + * @param x The point to check. + * @return True if the point is inside, false otherwise. + */ + bool EGS_DynamicGeometry::isInside(const EGS_Vector &x) { + EGS_Vector xt(x); + T.inverseTransform(xt); + return g->isInside(xt); + } + + /** + * @brief Determine the region of a point + * + * @param x The point to check. + * @return The region number (see EGS_BaseGeometry::isWhere() for details). + */ + int EGS_DynamicGeometry::isWhere(const EGS_Vector &x) { + EGS_Vector xt(x); + T.inverseTransform(xt); + return g->isWhere(xt); + } + + /** + * @brief Alias for EGS_DynamicGeometry::isWhere(). + * + * @param x The point to check. + * @return The region number (see EGS_BaseGeometry::isWhere() for details). + */ + int EGS_DynamicGeometry::inside(const EGS_Vector &x) { + return isWhere(x); + } + + /** + * @brief Get the medium index for a given region in the dynamic geometry. + * + * @param ireg The region index. + * @return The medium index. + */ + int EGS_DynamicGeometry::medium(int ireg) const { + return g->medium(ireg); + } + + /** + * @brief Calculate the distance to the outside of the dynamic geometry from a point along a direction. + * + * @param ireg The current region index. + * @param x The starting point. + * @param u The direction vector. + * @return The distance to the outside or 0 if the current region is invalid. + */ + EGS_Float EGS_DynamicGeometry::howfarToOutside(int ireg, const EGS_Vector &x, + const EGS_Vector &u) { + return ireg >= 0 ? g->howfarToOutside(ireg, x * T, u * T.getRotation()) : 0; + } + + /** + * @brief Calculate the distance to the next boundary in the dynamic geometry. + * + * @param ireg The current region index. + * @param x The starting point. + * @param u The direction vector. + * @param t The distance to the next boundary. + * @param newmed The medium index after boundary crossing. + * @param normal The normal vector at the boundary. + * @return The new region index after boundary crossing. + */ + int EGS_DynamicGeometry::howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u, + EGS_Float &t, int *newmed, EGS_Vector *normal) { + EGS_Vector xt(x), ut(u); + T.inverseTransform(xt); + T.rotateInverse(ut); + int inew = g->howfar(ireg, xt, ut, t, newmed, normal); + if (inew != ireg && normal) { + *normal = T.getRotation() * (*normal); + } + return inew; + } + + /** + * @brief Calculate the distance to the nearest boundary of the dynamic geometry from a point. + * + * @param ireg The current region index. + * @param x The point. + * @return The distance to the nearest boundary. + */ + EGS_Float EGS_DynamicGeometry::hownear(int ireg, const EGS_Vector &x) { + EGS_Vector xt(x); + T.inverseTransform(xt); + return g->hownear(ireg, xt); + } + + /** + * @brief Get the maximum step size allowed in the dynamic geometry. + * + * @return The maximum step size. + */ + int EGS_DynamicGeometry::getMaxStep() const { + return g->getMaxStep(); + } + + /** + * @brief Check if the specified region has a boolean property. + * + * @param ireg The region index. + * @param prop The boolean property to check. + * @return True if the region has the specified boolean property, false otherwise. + */ + bool EGS_DynamicGeometry::hasBooleanProperty(int ireg, EGS_BPType prop) const { + return g->hasBooleanProperty(ireg, prop); + } + + /** + * @brief Set a boolean property for the dynamic geometry. + * + * @param prop The boolean property to set. + */ + void EGS_DynamicGeometry::setBooleanProperty(EGS_BPType prop) { + g->setBooleanProperty(prop); + } + + /** + * @brief Add a boolean property to the dynamic geometry. + * + * @param bit The bit index of the boolean property to add. + */ + void EGS_DynamicGeometry::addBooleanProperty(int bit) { + g->addBooleanProperty(bit); + } + + /** + * @brief Set a boolean property for a range of regions in the dynamic geometry. + * + * @param prop The boolean property to set. + * @param start The starting region index. + * @param end The ending region index. + * @param step The step size between regions. + */ + void EGS_DynamicGeometry::setBooleanProperty(EGS_BPType prop, int start, int end, int step) { + g->setBooleanProperty(prop, start, end, step); + } + + /** + * @brief Add a boolean property to a range of regions in the dynamic geometry. + * + * @param bit The bit index of the boolean property to add. + * @param start The starting region index. + * @param end The ending region index. + * @param step The step size between regions. + */ + void EGS_DynamicGeometry::addBooleanProperty(int bit, int start, int end, int step) { + g->addBooleanProperty(bit, start, end, step); + } + + /** + * @brief Get the type of the dynamic geometry. + * + * @return The type as a string reference. + */ + const string &EGS_DynamicGeometry::getType() const { + return type; + } + + /** + * @brief Get the relative density scaling factor for a given region in the dynamic geometry. + * + * @param ireg The region index. + * @return The relative density scaling factor. + */ + EGS_Float EGS_DynamicGeometry::getRelativeRho(int ireg) const { + return g->getRelativeRho(ireg); + } + + /** + * @brief Get the boundary scaling factor for a given region in the dynamic geometry. + * + * @param ireg The region index. + * @return The boundary scaling factor. + */ + EGS_Float EGS_DynamicGeometry::getBScaling(int ireg) const { + return g->getBScaling(ireg); + } + + /** + * @brief Determine the next state of the dynamic geometry. + * + * This method obtains a time index, either from the simulation source or by + * sampling itself if the source yields nothing. It then obtains the + * corresponding position and orientation coordinates through the + * `getCoordGeom` method and creates and sets the dynamic geometry's + * transform. + * + * @param rndm Random number generator. + */ + void EGS_DynamicGeometry::getNextGeom(EGS_RandomGenerator *rndm) { + int errg = 1; + EGS_ControlPoint gipt; + + // Get the source from active application to extract time + EGS_Application *app = EGS_Application::activeApplication(); + while (errg) { + // Get time from source if it exists (otherwise gives -1) + ptime = app->getTimeIndex(); + if (ptime < 0) { + // If no time is given by the source, randomly sample from 0 to 1 + ptime = rndm->getUniform(); + // Set randomly sampled time index for all objects in the + // simulation (through base source) + app->setTimeIndex(ptime); + } + // Run the `getCoordGeom` method that will sample the control points + // given to find the transformation that will be applied for the + // current history + errg = getCoordGeom(ptime, gipt); + } + + // Create and set the current geometry transformation using the sampled + // coordinates from `getCoordGeom`. This is where the overloaded + // `EGS_AffineTransform` is used + EGS_AffineTransform *tDG = EGS_AffineTransform::getTransformation(gipt.trnsl, gipt.rot); + setTransformation(*tDG); + + // Call `getNextGeom` on base geometry in case there are lower-level + // dynamic geometries + g->getNextGeom(rndm); + } + + /** + * @brief Update the position of the dynamic geometry. + * + * This method is used to update the geometry state as needed in egs_view. + * It takes the desired time index, computes the corresponding coordinates, + * and sets the geometry transform to update the egs_view display. + * + * @param time Desired time index for the update. + */ + void EGS_DynamicGeometry::updatePosition(EGS_Float time) { + int errg = 1; + EGS_ControlPoint gipt; + + // Run the `getCoordGeom` method that will use the control points given + // to find the transformation that will be applied for the current time + // index + errg = getCoordGeom(time, gipt); + + // Create and set the geometry transform with the updated coordinates + EGS_AffineTransform *tDG = EGS_AffineTransform::getTransformation(gipt.trnsl, gipt.rot); + setTransformation(*tDG); + + // Call `updatePosition` on the base to allow lower-level geometries to + // update as needed + g->updatePosition(time); + } + + /** + * @brief Check if the simulation geometry contains dynamic geometry. + * + * This method is used to determine whether the simulation geometry contains + * dynamic geometry. It sets the `hasdynamic` flag to true if the simulation + * contains dynamic geometry. + * + * @param hasdynamic Boolean flag to indicate if dynamic geometry is + * present. + */ + void EGS_DynamicGeometry::containsDynamic(bool &hasdynamic) { + // If the dynamic geometry implementation of `containsDynamic` is + // called, the simulation does indeed contain a dynamic geometry, so the + // boolean flag is set to true + hasdynamic = true; + } + +} diff --git a/HEN_HOUSE/egs++/geometry/egs_dynamic_geometry/egs_dynamic_geometry.h b/HEN_HOUSE/egs++/geometry/egs_dynamic_geometry/egs_dynamic_geometry.h new file mode 100644 index 000000000..9f1688df6 --- /dev/null +++ b/HEN_HOUSE/egs++/geometry/egs_dynamic_geometry/egs_dynamic_geometry.h @@ -0,0 +1,438 @@ +/* +############################################################################### +# +# EGSnrc egs++ dynamic_geometry geometry 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: Alexandre Demelo 2023 +# +# Contributors: Reid Townson +# +############################################################################### +*/ + +/*! \file egs_dynamic_geometry.h + * \brief A dynamic geometry: header + * \RT + */ + +#ifndef EGS_DYNAMIC_GEOMETRY_ +#define EGS_DYNAMIC_GEOMETRY_ + +#include "egs_base_geometry.h" +#include "egs_transformations.h" +#include "egs_rndm.h" +#include "egs_vector.h" +#include "egs_application.h" +#include +#include +#include + +#ifdef WIN32 + + #ifdef BUILD_DYNAMIC_GEOMETRY_DLL + #define EGS_DYNAMIC_GEOMETRY_EXPORT __declspec(dllexport) + #else + #define EGS_DYNAMIC_GEOMETRY_EXPORT __declspec(dllimport) + #endif + #define EGS_DYNAMIC_GEOMETRY_LOCAL + +#else + + #ifdef HAVE_VISIBILITY + #define EGS_DYNAMIC_GEOMETRY_EXPORT __attribute__ ((visibility ("default"))) + #define EGS_DYNAMIC_GEOMETRY_LOCAL __attribute__ ((visibility ("hidden"))) + #else + #define EGS_DYNAMIC_GEOMETRY_EXPORT + #define EGS_DYNAMIC_GEOMETRY_LOCAL + #endif + +#endif + +/*! \brief A dynamic geometry. + +\ingroup Geometry \ingroup CompositeG + +An dynamic geometry is a geometry that takes a random point from another +geometry and then applies a transformation, using a time sampling and +interpolating between control points. + +An dynamic geometry is defined using \verbatim :start geometry: name = +... library = egs_dynamic_geometry my geometry = name of a predefined +geometry that we want to add motion to :start motion: # units of cm and degrees +control point = time(1) xtrans(1) ytrans(1) ztrans(1) xrot(1) yrot(1) zrot(1) +control point = time(2) xtrans(2) ytrans(2) ztrans(2) xrot(2) yrot(2) zrot(2) + . + . + . + control point = time(N) xtrans(N) ytrans(N) ztrans(N) xrot(N) yrot(N) zrot(N) + :stop motion: +:stop geometry: \endverbatim + +Control points must be defined such that time(i+1)>=time(i), where time(i) is +the value of time for control point i. The time(i) are automatically normalized +by time(N), where N is the number of control points. + +A translation from the starting position of the geometry is applied according to +x, y and z. A rotation follows the same rotation technique as in +EGS_AffineTransform, using the rotation input parameter for 2 or 3 values. +Angles are in degrees and translations in cm. + +Continuous, dynamic motion between control points is simulated by choosing a +random number, R, on (0,1] and, for time(i)0.0, in the case where a user desires to eliminate particles associated +with a range of time values, but there will be a lot of warning messages. + +*/ + +class EGS_DYNAMIC_GEOMETRY_EXPORT EGS_DynamicGeometry : + public EGS_BaseGeometry { + +public: + + /*! + * \struct EGS_ControlPoint + * \brief Structure to store control point information for dynamic geometry. + */ + struct EGS_ControlPoint { + EGS_Float time; //!< Time index for control point + vector trnsl; //!< Vector specifying x, y, z translation + vector rot; //!< Vector specifying rotation + }; + + /*! \brief Construct a dynamic geometry using \a G as the geometry and \a cpts as the control points. + * + * \param G Pointer to the base geometry to be transformed dynamically. + * \param dyninp Input containing dynamic geometry specifications. + * \param Name Name of the dynamic geometry. + */ + EGS_DynamicGeometry(EGS_BaseGeometry *G, EGS_Input *dyninp, const string &Name = "") : EGS_BaseGeometry(Name), g(G) { + type = g->getType() + "D"; + nreg = g->regions(); + is_convex = g->isConvex(); + has_rho_scaling = g->hasRhoScaling(); + has_B_scaling = g->hasBScaling(); + EGS_DynamicGeometry::buildDynamicGeometry(g, dyninp); + + if (cpts.size() < 2) { + egsWarning("EGS_DynamicGeometry: not enough or missing control points.\n"); + } + else { + if (cpts[0].time > 0.0) { + egsWarning("EGS_DynamicGeometry: time index of control point 1 > 0.0. This will generate many warning messages.\n"); + } + int npts = cpts.size(); + for (int i = 0; i < npts; i++) { + if (i > 0 && cpts[i].time < cpts[i - 1].time - epsilon) { + egsWarning("EGS_DynamicGeometry: time index of control point %i < time index of control point %i\n", i, i - 1); + } + if (cpts[i].time < 0.0) { + egsWarning("EGS_DynamicGeometry: time index of control point %i < 0.0\n", i); + } + } + // Normalize time values + for (int i = 0; i < npts - 1; i++) { + cpts[i].time /= cpts[npts - 1].time; + } + } + + }; + + /*! Destructor. */ + ~EGS_DynamicGeometry() { + if (!g->deref()) { + delete g; + } + }; + + /*! + * Sets the current state transform of the geometry. This is called when + * checking location. Same as gtransformed. + * + * \param t Affine transformation to set. + */ + void setTransformation(EGS_AffineTransform t) { + T = t; + }; + + /*! + * Computes intersections of a particle with the dynamic geometry. + * + * \param ireg Region index. + * \param n Number of particles. + * \param x Particle position. + * \param u Particle direction. + * \param isections Geometry intersections. + * \return Number of intersections. + */ + int computeIntersections(int ireg, int n, const EGS_Vector &x, + const EGS_Vector &u, EGS_GeometryIntersections *isections); + + /*! + * Checks if a region is real. + * + * \param ireg Region index. + * \return True if the region is real, false otherwise. + */ + bool isRealRegion(int ireg) const; + + /*! + * Checks if a point is inside the dynamic geometry. + * + * \param x Point to check. + * \return True if the point is inside, false otherwise. + */ + bool isInside(const EGS_Vector &x); + + /*! + * Checks the location of a point. + * + * \param x Point to check. + * \return Index of the region containing the point. + */ + int isWhere(const EGS_Vector &x); + + /*! + * Alias for isWhere method. + * + * \param x Point to check. + * \return Index of the region containing the point. + */ + int inside(const EGS_Vector &x); + + /*! + * Returns the medium index of a region. + * + * \param ireg Region index. + * \return Medium index. + */ + int medium(int ireg) const; + + /*! + * Computes the distance to the outside of the dynamic geometry. + * + * \param ireg Region index. + * \param x Particle position. + * \param u Particle direction. + * \return Distance to the outside. + */ + EGS_Float howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u); + + /*! + * Computes the distance to the nearest boundary. + * + * \param ireg Region index. + * \param x Particle position. + * \param u Particle direction. + * \param t Output: distance to the nearest boundary. + * \param newmed Output: new medium index after crossing the boundary. + * \param normal Output: normal vector at the intersection point. + * \return Index of the new region. + */ + int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u, EGS_Float &t, int *newmed=0, EGS_Vector *normal=0); + + /*! + * Computes the distance to the nearest boundary. + * + * \param ireg Region index. + * \param x Particle position. + * \return Distance to the nearest boundary. + */ + EGS_Float hownear(int ireg, const EGS_Vector &x); + + /*! + * Returns the maximum step allowed for the dynamic geometry. + * + * \return Maximum allowed step. + */ + int getMaxStep() const; + + /*! + * Checks if the dynamic geometry has a specific boolean property in a region. + * + * \param ireg Region index. + * \param prop Boolean property type. + * \return True if the region has the specified property, false otherwise. + */ + bool hasBooleanProperty(int ireg, EGS_BPType prop) const; + + /*! + * Sets a boolean property for the dynamic geometry. + * + * \param prop Boolean property type. + */ + void setBooleanProperty(EGS_BPType prop); + + /*! + * Adds a boolean property for the dynamic geometry. + * + * \param bit Bit index of the property to add. + */ + void addBooleanProperty(int bit); + + /*! + * Sets a boolean property for a range of regions in the dynamic geometry. + * + * \param prop Boolean property type. + * \param start Start index of the range. + * \param end End index of the range. + * \param step Step size for the range. + */ + void setBooleanProperty(EGS_BPType prop, int start, int end, int step=1); + + /*! + * Adds a boolean property for a range of regions in the dynamic geometry. + * + * \param bit Bit index of the property to add. + * \param start Start index of the range. + * \param end End index of the range. + * \param step Step size for the range. + */ + void addBooleanProperty(int bit, int start, int end, int step=1); + + /*! + * Returns the type of the dynamic geometry. + * + * \return Type of the dynamic geometry. + */ + const string &getType() const; + + /*! + * Gets the relative density of a region in the dynamic geometry. + * + * \param ireg Region index. + * \return Relative density of the region. + */ + EGS_Float getRelativeRho(int ireg) const; + + /*! + * Sets the relative density for a range of regions in the dynamic geometry. + * + * \param start Start index of the range. + * \param end End index of the range. + * \param rho Relative density to set. + */ + void setRelativeRho(int start, int end, EGS_Float rho); + + /*! + * Sets the relative density for a range of regions in the dynamic geometry + * using input specifications. + * + * \param input Input containing relative density specifications. + */ + void setRelativeRho(EGS_Input *); + + /*! + * Gets the magnetic field scaling factor for a region in the dynamic + * geometry. + * + * \param ireg Region index. \return Magnetic field scaling factor. + */ + EGS_Float getBScaling(int ireg) const; + + /*! + * Sets the magnetic field scaling factor for a range of regions in the dynamic geometry. + * + * \param start Start index of the range. + * \param end End index of the range. + * \param bf Magnetic field scaling factor to set. + */ + void setBScaling(int start, int end, EGS_Float bf); + + /*! + * Sets the magnetic field scaling factor for a range of regions in the + * dynamic geometry using input specifications. + * + * \param input Input containing magnetic field scaling factor + * specifications. + */ + void setBScaling(EGS_Input *); + + /*! + * Retrieves regions labeled with a given string. + * + * \param str Label to search for. + * \param regs Output: List of region indices with the specified label. + */ + void getLabelRegions(const string &str, vector ®s); + + /*! + * Updates the next particle state for geometries. It is tasked with + * determining the next state of the dynamic geometry. + * + * \param rndm Random number generator. + */ + void getNextGeom(EGS_RandomGenerator *rndm); + + /*! + * Updates the position of the dynamic geometry to the specified time. + * + * \param time Time index to update to. + */ + void updatePosition(EGS_Float time); + + /*! + * Determines whether the simulation geometry contains a dynamic geometry. + * + * \param hasdynamic Output: True if the simulation contains a dynamic geometry, false otherwise. + */ + void containsDynamic(bool &hasdynamic); + +protected: + EGS_BaseGeometry *g; //!< The geometry undergoing dynamic motion + string type; //!< The geometry type + EGS_AffineTransform T; //!< Affine transformation representing the current state + vector cpts; //!< Control points for dynamic motion + int ncpts; //!< Number of control points + EGS_Float ptime; //!< Time index corresponding to the particle + + /*! + * \brief Don't define media in the transformed geometry definition. + * + * This function is re-implemented to warn the user not to define media in + * the definition of a transformed geometry. Instead, media should be + * defined when specifying the geometry to be transformed. + */ + void setMedia(EGS_Input *inp, int, const int *); + + /*! + * \brief Extract coordinates for the next dynamic geometry position. + * + * \param rand Random number for time sampling. + * \param gipt EGS_ControlPoint structure to store the coordinates. + * \return 0 if successful, otherwise 1. + */ + int getCoordGeom(EGS_Float rand, EGS_ControlPoint &gipt); + + /*! + * \brief Builds the dynamic geometry using input specifications. + * + * \param g Pointer to the base geometry to be transformed dynamically. + * \param dyninp Input containing dynamic geometry specifications. + */ + void buildDynamicGeometry(EGS_BaseGeometry *g, EGS_Input *dyninp); +}; + +#endif diff --git a/HEN_HOUSE/egs++/geometry/egs_genvelope/egs_envelope_geometry.cpp b/HEN_HOUSE/egs++/geometry/egs_genvelope/egs_envelope_geometry.cpp index f94a7f468..02d7dc917 100644 --- a/HEN_HOUSE/egs++/geometry/egs_genvelope/egs_envelope_geometry.cpp +++ b/HEN_HOUSE/egs++/geometry/egs_genvelope/egs_envelope_geometry.cpp @@ -319,7 +319,7 @@ void EGS_EnvelopeGeometry::printInfo() const { g->getType().c_str()); egsInformation(" inscribed geometries:\n"); for (int j=0; jgetName().c_str(),geometries[j]->getType().c_str()); + geometries[j]->getName().c_str(),geometries[j]->getType().c_str()); egsInformation( "=======================================================\n"); } @@ -330,7 +330,7 @@ void EGS_FastEnvelope::printInfo() const { g->getType().c_str()); egsInformation(" inscribed geometries:\n"); for (int j=0; jgetName().c_str(),geometries[j]->getType().c_str()); + geometries[j]->getName().c_str(),geometries[j]->getType().c_str()); egsInformation( "=======================================================\n"); } diff --git a/HEN_HOUSE/egs++/geometry/egs_genvelope/egs_envelope_geometry.h b/HEN_HOUSE/egs++/geometry/egs_genvelope/egs_envelope_geometry.h index 184258883..6ee53ba4d 100644 --- a/HEN_HOUSE/egs++/geometry/egs_genvelope/egs_envelope_geometry.h +++ b/HEN_HOUSE/egs++/geometry/egs_genvelope/egs_envelope_geometry.h @@ -27,6 +27,7 @@ # Frederic Tessier # Reid Townson # Max Orok +# Alexandre Demelo # ############################################################################### */ @@ -523,6 +524,40 @@ class EGS_ENVELOPEG_EXPORT EGS_EnvelopeGeometry : public EGS_BaseGeometry { return g->hownear(ireg,x); }; + + void getNextGeom(EGS_RandomGenerator *rndm) { + // calls getNextGeom on its component geometries to update dynamic + // geometries in the simulation + for (int j=0; jgetNextGeom(rndm); + } + g->getNextGeom(rndm); + }; + + void updatePosition(EGS_Float time) { + // calls updatePosition on its component geometries to update dynamic + // geometries in the simulation + for (int j=0; jupdatePosition(time); + } + g->updatePosition(time); + }; + + void containsDynamic(bool &hasdynamic) { + // calls containsDynamic on its component geometries (only calls if + // hasDynamic is false, as if it is true we already found one) + for (int j=0; jcontainsDynamic(hasdynamic); + } + } + if (!hasdynamic) { + g->containsDynamic(hasdynamic); + } + }; + + + int getMaxStep() const { int nstep = g->getMaxStep(); for (int j=0; jgetNextGeom(rndm); + } + g->getNextGeom(rndm); + }; + + void updatePosition(EGS_Float time) { + // calls updatePosition on its component geometries to update dynamic + // geometries in the simulation + for (int j=0; jupdatePosition(time); + } + g->updatePosition(time); + }; + + void containsDynamic(bool &hasdynamic) { + // calls containsDynamic on its component geometries (only calls if + // hasDynamic is false, as if it is true we already found one) + for (int j=0; jcontainsDynamic(hasdynamic); + } + } + if (!hasdynamic) { + g->containsDynamic(hasdynamic); + } + }; + + int inside(const EGS_Vector &x) { return isWhere(x); }; diff --git a/HEN_HOUSE/egs++/geometry/egs_gtransformed/egs_gtransformed.h b/HEN_HOUSE/egs++/geometry/egs_gtransformed/egs_gtransformed.h index 83247d9ee..4c7d0d8ea 100644 --- a/HEN_HOUSE/egs++/geometry/egs_gtransformed/egs_gtransformed.h +++ b/HEN_HOUSE/egs++/geometry/egs_gtransformed/egs_gtransformed.h @@ -26,6 +26,7 @@ # Contributors: Frederic Tessier # Reid Townson # Ernesto Mainegra-Hing +# Alexandre Demelo # ############################################################################### */ @@ -217,6 +218,27 @@ class EGS_GTRANSFORMED_EXPORT EGS_TransformedGeometry : //return g->hownear(ireg,x*T); }; + void getNextGeom(EGS_RandomGenerator *rndm) { + // calls getNextGeom on its component geometries to update dynamic + // geometries in the simulation + g->getNextGeom(rndm); + + }; + + void updatePosition(EGS_Float time) { + // calls updatePosition on its component geometries to update dynamic + // geometries in the simulation + g->updatePosition(time); + }; + + void containsDynamic(bool &hasdynamic) { + // calls containsDynamic on its component geometries (only calls if + // hasDynamic is false, as if it is true we already found one) + if (!hasdynamic) { + g->containsDynamic(hasdynamic); + } + }; + int getMaxStep() const { return g->getMaxStep(); }; diff --git a/HEN_HOUSE/egs++/geometry/egs_mesh/egs_mesh.h b/HEN_HOUSE/egs++/geometry/egs_mesh/egs_mesh.h index 9491d0583..e8b56565e 100644 --- a/HEN_HOUSE/egs++/geometry/egs_mesh/egs_mesh.h +++ b/HEN_HOUSE/egs++/geometry/egs_mesh/egs_mesh.h @@ -365,9 +365,9 @@ class EGS_MESH_EXPORT EGS_Mesh : public EGS_BaseGeometry { const auto &node_indices = elt_node_indices_.at(element); return Nodes { nodes_.at(node_indices[0]), - nodes_.at(node_indices[1]), - nodes_.at(node_indices[2]), - nodes_.at(node_indices[3]) + nodes_.at(node_indices[1]), + nodes_.at(node_indices[2]), + nodes_.at(node_indices[3]) }; } diff --git a/HEN_HOUSE/egs++/geometry/egs_mesh/mesh_neighbours.h b/HEN_HOUSE/egs++/geometry/egs_mesh/mesh_neighbours.h index bfbb76c77..da42e4610 100644 --- a/HEN_HOUSE/egs++/geometry/egs_mesh/mesh_neighbours.h +++ b/HEN_HOUSE/egs++/geometry/egs_mesh/mesh_neighbours.h @@ -159,8 +159,8 @@ SharedNodes elements_around_nodes(const std::vector> tetrahedron_neighbours( - const std::vector &elements, -egs_mesh::internal::PercentCounter &progress) { + const std::vector &elements, + egs_mesh::internal::PercentCounter &progress) { progress.start(elements.size()); const std::size_t NUM_FACES = 4; const auto shared_nodes = mesh_neighbours::internal::elements_around_nodes(elements); diff --git a/HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.cpp b/HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.cpp index 8daf2ca5a..330ab00bd 100644 --- a/HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.cpp +++ b/HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.cpp @@ -1342,7 +1342,7 @@ extern "C" { "construct a ND geometry with a single dimension?\n"); input->print(0,cerr); for (int j=0; jisConvex()) { diff --git a/HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.h b/HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.h index 15ac5c716..ce1f5b6ce 100644 --- a/HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.h +++ b/HEN_HOUSE/egs++/geometry/egs_nd_geometry/egs_nd_geometry.h @@ -28,6 +28,7 @@ # Frederic Tessier # Reid Townson # Randle Taylor +# Alexandre Demelo # ############################################################################### */ @@ -554,6 +555,32 @@ class EGS_NDG_EXPORT EGS_NDGeometry : public EGS_BaseGeometry { return nstep; }; + void getNextGeom(EGS_RandomGenerator *rndm) { + // calls getNextGeom on its component geometries to update dynamic + // geometries in the simulation + for (int j=0; jgetNextGeom(rndm); + } + }; + + void updatePosition(EGS_Float time) { + // calls updatePosition on its component geometries to update dynamic + // geometries in the simulation + for (int j=0; jupdatePosition(time); + } + }; + + void containsDynamic(bool &hasdynamic) { + // calls containsDynamic on its component geometries (only calls if + // hasDynamic is false, as if it is true we already found one) + for (int j=0; jcontainsDynamic(hasdynamic); + } + } + }; + const string &getType() const { return type; }; diff --git a/HEN_HOUSE/egs++/geometry/egs_planes/egs_planes.cpp b/HEN_HOUSE/egs++/geometry/egs_planes/egs_planes.cpp index 24e4ed8d0..962f1fbce 100644 --- a/HEN_HOUSE/egs++/geometry/egs_planes/egs_planes.cpp +++ b/HEN_HOUSE/egs++/geometry/egs_planes/egs_planes.cpp @@ -129,11 +129,11 @@ extern "C" { } EGS_BaseGeometry *g; if (type == "EGS_Xplanes") g = new EGS_PlanesX(pos,"", - EGS_XProjector(xproj_type)); + EGS_XProjector(xproj_type)); else if (type == "EGS_Yplanes") g = new EGS_PlanesY(pos,"", - EGS_YProjector(yproj_type)); + EGS_YProjector(yproj_type)); else if (type == "EGS_Zplanes") g = new EGS_PlanesZ(pos,"", - EGS_ZProjector(zproj_type)); + EGS_ZProjector(zproj_type)); else if (type == "EGS_Planes") { vector a; err = input->getInput("normal",a); diff --git a/HEN_HOUSE/egs++/geometry/egs_roundrect_cylinders/egs_roundrect_cylinders.cpp b/HEN_HOUSE/egs++/geometry/egs_roundrect_cylinders/egs_roundrect_cylinders.cpp index 626f5437f..3394f2074 100644 --- a/HEN_HOUSE/egs++/geometry/egs_roundrect_cylinders/egs_roundrect_cylinders.cpp +++ b/HEN_HOUSE/egs++/geometry/egs_roundrect_cylinders/egs_roundrect_cylinders.cpp @@ -137,13 +137,13 @@ extern "C" { EGS_BaseGeometry *g; if (type == "EGS_RoundRectCylindersXY") g = new EGS_RoundRectCylindersT(x_widths,y_widths,radii,xo,"", - EGS_XProjector("X"),EGS_YProjector("Y")); + EGS_XProjector("X"),EGS_YProjector("Y")); else if (type == "EGS_RoundRectCylindersXZ") g = new EGS_RoundRectCylindersT(x_widths,y_widths,radii,xo,"", - EGS_XProjector("X"),EGS_ZProjector("Z")); + EGS_XProjector("X"),EGS_ZProjector("Z")); else if (type == "EGS_RoundRectCylindersYZ") g = new EGS_RoundRectCylindersT(x_widths,y_widths,radii,xo,"", - EGS_YProjector("Y"),EGS_ZProjector("Z")); + EGS_YProjector("Y"),EGS_ZProjector("Z")); else { vector ax, ay; err = input->getInput("x-axis",ax); @@ -159,8 +159,8 @@ extern "C" { return 0; } g = new EGS_RoundRectCylindersT(x_widths,y_widths,radii,xo,"", - EGS_Projector(EGS_Vector(ax[0],ax[1],ax[2]),"Any"), - EGS_Projector(EGS_Vector(ay[0],ay[1],ay[2]),"")); + EGS_Projector(EGS_Vector(ax[0],ax[1],ax[2]),"Any"), + EGS_Projector(EGS_Vector(ay[0],ay[1],ay[2]),"")); } g->setName(input); g->setLabels(input); diff --git a/HEN_HOUSE/egs++/geometry/egs_smart_envelope/egs_smart_envelope.cpp b/HEN_HOUSE/egs++/geometry/egs_smart_envelope/egs_smart_envelope.cpp index b0e72b53f..c40416f79 100644 --- a/HEN_HOUSE/egs++/geometry/egs_smart_envelope/egs_smart_envelope.cpp +++ b/HEN_HOUSE/egs++/geometry/egs_smart_envelope/egs_smart_envelope.cpp @@ -223,9 +223,9 @@ void EGS_SmartEnvelope::printInfo() const { g->getType().c_str()); egsInformation(" inscribed geometries:\n"); for (int j=0; jgetName().c_str(),geometries[j]->getType().c_str(), - reg_to_base[j],(int)itype[j]); + " itype=%d\n", + geometries[j]->getName().c_str(),geometries[j]->getType().c_str(), + reg_to_base[j],(int)itype[j]); egsInformation( "=======================================================\n"); } diff --git a/HEN_HOUSE/egs++/geometry/egs_smart_envelope/egs_smart_envelope.h b/HEN_HOUSE/egs++/geometry/egs_smart_envelope/egs_smart_envelope.h index 912a81ba7..b449712b7 100644 --- a/HEN_HOUSE/egs++/geometry/egs_smart_envelope/egs_smart_envelope.h +++ b/HEN_HOUSE/egs++/geometry/egs_smart_envelope/egs_smart_envelope.h @@ -25,6 +25,7 @@ # # Contributors: Frederic Tessier # Ernesto Mainegra-Hing +# Alexandre Demelo # ############################################################################### # @@ -317,6 +318,37 @@ class EGS_SMART_ENVELOPE_EXPORT EGS_SmartEnvelope : public EGS_BaseGeometry { return g->hownear(ireg,x); }; + void getNextGeom(EGS_RandomGenerator *rndm) { + // calls getNextGeom on its component geometries to update dynamic + // geometries in the simulation + for (int j=0; jgetNextGeom(rndm); + } + g->getNextGeom(rndm); + }; + + void updatePosition(EGS_Float time) { + // calls updatePosition on its component geometries to update dynamic + // geometries in the simulation + for (int j=0; jupdatePosition(time); + } + g->updatePosition(time); + }; + + void containsDynamic(bool &hasdynamic) { + // calls containsDynamic on its component geometries (only calls if + // hasDynamic is false, as if it is true we already found one) + for (int j=0; jcontainsDynamic(hasdynamic); + } + } + if (!hasdynamic) { + g->containsDynamic(hasdynamic); + } + }; + int getMaxStep() const { int nstep = g->getMaxStep() + n_in; for (int j=0; jgetNextGeom(rndm); + } + }; + + void updatePosition(EGS_Float time) { + // calls updatePosition on its component geometries to update dynamic + // geometries in the simulation + for (int j=0; jupdatePosition(time); + } + }; + + void containsDynamic(bool &hasdynamic) { + // calls containsDynamic on its component geometries (only calls if + // hasDynamic is false, as if it is true we already found one) + for (int j=0; jcontainsDynamic(hasdynamic); + } + } + }; + int getMaxStep() const { int nstep = 1; for (int j=0; jsetName(input); result->setBoundaryTolerance(input); result->setLabels(input); diff --git a/HEN_HOUSE/egs++/shapes/egs_dynamic_shape/Makefile b/HEN_HOUSE/egs++/shapes/egs_dynamic_shape/Makefile new file mode 100644 index 000000000..533705833 --- /dev/null +++ b/HEN_HOUSE/egs++/shapes/egs_dynamic_shape/Makefile @@ -0,0 +1,44 @@ + +############################################################################### +# +# EGSnrc egs++ makefile to build dynamic shape +# 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: Reid Townson +# +# Contributors: +# +############################################################################### + + +include $(EGS_CONFIG) +include $(SPEC_DIR)egspp.spec +include $(SPEC_DIR)egspp_$(my_machine).conf + +DEFS = $(DEF1) -DBUILD_DYNAMIC_SHAPE_DLL + +library = egs_dynamic_shape +lib_files = egs_dynamic_shape +my_deps = $(common_shape_deps) +extra_dep = $(addprefix $(DSOLIBS), $(my_deps)) + +include $(SPEC_DIR)egspp_libs.spec + +$(make_depend) diff --git a/HEN_HOUSE/egs++/shapes/egs_dynamic_shape/egs_dynamic_shape.cpp b/HEN_HOUSE/egs++/shapes/egs_dynamic_shape/egs_dynamic_shape.cpp new file mode 100644 index 000000000..ae99670fe --- /dev/null +++ b/HEN_HOUSE/egs++/shapes/egs_dynamic_shape/egs_dynamic_shape.cpp @@ -0,0 +1,289 @@ +/* +############################################################################### +# +# EGSnrc egs++ dynamic shape +# 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: Reid Townson +# +# Contributors: +# +############################################################################### +*/ + + +/*! \file egs_dynamic_shape.cpp + * \brief Implementation of a dynamic shape + * \RT + */ + +#include "egs_dynamic_shape.h" +#include "egs_input.h" +#include "egs_functions.h" +#include + +extern "C" { + + /*! + * \brief Factory function to create an instance of EGS_DynamicShape + * \param input Input specifying the dynamic shape + * \param f EGS_ObjectFactory pointer + * \return Pointer to the created EGS_DynamicShape instance + */ + EGS_DYNAMIC_SHAPE_EXPORT EGS_BaseShape *createShape(EGS_Input *input, + EGS_ObjectFactory *f) { + if (!input) { + egsWarning("createShape(dynamic shape): null input?\n"); + return 0; + } + EGS_Input *ishape = input->takeInputItem("shape",false); + EGS_BaseShape *shape=0;; + if (ishape) { + shape = EGS_BaseShape::createShape(ishape); + delete ishape; + } + if (!shape) { + string shape_name; + int err = input->getInput("shape name",shape_name); + if (err) { + egsWarning("createShape(dynamic shape): no inline shape definition" + " and no 'shape name' keyword\n"); + return 0; + } + shape = EGS_BaseShape::getShape(shape_name); + if (!shape) { + egsWarning("createShape(dynamic shape): no shape named %s " + "exists\n",shape_name.c_str()); + return 0; + } + } + + // Now getting motion information for dynamic component + EGS_Input *dyninp = input->takeInputItem("motion"); + + if (dyninp) { + EGS_DynamicShape *s = new EGS_DynamicShape(shape, dyninp, "", f); + s->setName(input); + return s; + } + else { + egsWarning("EGS_DynamicShape: no control points input.\n"); + return 0; + } + } + + /*! + * \brief Builds the dynamic shape using input specifications + * \param dyninp Input containing dynamic shape specifications + */ + void EGS_DynamicShape::buildDynamicShape(EGS_Input *dyninp) { + stringstream itos; + ncpts=0; + vector point; + EGS_ControlPoint cpt; + int err; + int icpts=1; + itos << icpts; + string inputTag = "control point"; + string inputTag_backCompat = "control point " + itos.str(); + EGS_Input *currentInput; + int rotsize=0; + + // Control points read one by one from motion block in dynamic geometry definition, and saved as a vector to points + while (true) { + currentInput = dyninp->takeInputItem(inputTag); + if (!currentInput || currentInput->getInput(inputTag, point)) { + currentInput = dyninp->takeInputItem(inputTag_backCompat); + if (!currentInput || currentInput->getInput(inputTag_backCompat, point)) { + delete currentInput; + break; + } + } + delete currentInput; + + // Checking the size to make sure it is a valid control point input + if (point.size()!=6 && point.size()!=7) { + egsFatal("EGS_DynamicShape: Control point %i must specify either 6 or 7 values.\n",icpts); + } + else { + if (ncpts == 0) { + rotsize=point.size();// Variable to make sure all control point definitions have consistent formats + } + if (ncpts>0 && point[0] < cpts[ncpts-1].time) {// Make sure each time index is larger than the last + egsFatal("EGS_DynamicShape: Time index of control point %i < time index of control point %i\n",icpts,ncpts); + } + else if (point[0] < 0.) {// Checks that time index is valid (larger than zero) + egsFatal("EGS_DynamicShape: Time index of control point %i < 0.0\n",icpts); + } + else if (ncpts> 0 && point.size() != rotsize) { // Checks that control point formats follow the first + egsFatal("EGS_DynamicShape: Rotation definition inconsistent \n"); + } + else { + ncpts++; + + if (ncpts ==1 && point[0] > 0.0) { + egsWarning("EGS_DynamicShape: Time index of control point 1 > 0.0. This will generate many warning messages.\n"); + } + + vector T_vect;// Vector storing translation information + vector R_vect;// Vector storing rotation information + // Add translation coordinates to the translation vector + T_vect.push_back(point[1]); + T_vect.push_back(point[2]); + T_vect.push_back(point[3]); + // Add rotation coordinates to rotation vector (cases differentiate cpt formats, 6 is 2 rotation parameters, 7 is 3) + // In each case vector order is determined by format of EGS_rotationMatrix constructor (EGS_Float arguments) + if (point.size()==6) { + R_vect.push_back(point[6]);//rotation about z + R_vect.push_back(point[4]);//rotation about x + }//2 number case + if (point.size()==7) { + R_vect.push_back(point[4]);//rotation about x + R_vect.push_back(point[5]);//rotation about y + R_vect.push_back(point[6]);//rotation about z + }//3 number case + + // Add it to the vector of control points + cpt.time = point[0]; + cpt.trnsl = T_vect; + cpt.rot = R_vect; + cpts.push_back(cpt);// Add created control point to list of control points + icpts++; + itos.str(""); + itos << icpts; + inputTag_backCompat = "control point " + itos.str();// Define next control point i string for getInput in while condition + } + } + } + if (ncpts<=1) { + egsFatal("EGS_DynamicShape: not enough or missing control points.\n"); + } + if (cpts[ncpts-1].time == 0.0) { + egsFatal("EGS_DynamicShape: time index of last control point = 0. Something's wrong.\n"); + } + else { + // Normalize time index to max. value + for (int i=0; i<=ncpts-1; i++) {// I changed the normalization here (and in dynamic source) to have <= instead of <. Otherwise last cpt never gets normalized, which is not an issue for final cpt time>1 but is an issue if final cpt time<1. + cpts[i].time /= cpts[ncpts-1].time; + } + } + updatePosition(0); // Sets position to initial time in egs_view upon opening + }; + + /*! + * \brief Get the next state of the dynamic shape + * \param rndm Random number generator + */ + void EGS_DynamicShape::getNextShapePosition(EGS_RandomGenerator *rndm) { + int errg = 1; + EGS_ControlPoint gipt; + + // Here get source from activeapplication in order to extract time + EGS_Application *app = EGS_Application::activeApplication(); + while (errg) { + // Gets time if it's already set (otherwise gives -1). + ptime = app->getTimeIndex(); + + if (ptime<0) { + // If no time is given by the source the shape will randomly sample from 0 to 1. + ptime = rndm->getUniform(); + + // Set randomly sampled time index for all objects in the simulation + app->setTimeIndex(ptime); + } + + // Now run the get coord method that will sample the cpt given to find the transformation that will be applied for the current history + errg = getCoord(ptime,gipt); + } + + // Create and set the current shape transformation using the sampled coordinates from getCoord. This is where overloaded EGS_AffineTransform is used + EGS_AffineTransform *tDG = EGS_AffineTransform::getTransformation(gipt.trnsl, gipt.rot); + shape->setTransformation(tDG); + delete tDG; + + // Call getNextShapePosition on base shape in case there are lower level dynamic shapes + shape->getNextShapePosition(rndm); + }; + + /*! + * \brief Extract coordinates for the next dynamic shape position + * \param rand Random number for time sampling + * \param gipt EGS_ControlPoint structure to store the coordinates + * \return 0 if successful, otherwise 1 + */ + int EGS_DynamicShape::getCoord(EGS_Float rand, EGS_ControlPoint &gipt) { + int iindex=0; + int i; + + // The following loop determines which 2 control points the current time index falls between + for (i=0; i translation_LB=cpts[iindex-1].trnsl; + vector translation_UB=cpts[iindex].trnsl; + vector translation_samp; + + // The following is a factor (between 0 and 1) used for sampling. Essentially, it represents the fractional position within the interval + EGS_Float factor = (rand-cpts[iindex-1].time)/(cpts[iindex].time-cpts[iindex-1].time); + + // Translations are given as the lower bound plus the length of the interval multiplied by the fractional position factor. So its essentially lower bound + n% of the interval length + translation_samp.push_back(translation_LB[0]+(translation_UB[0]-translation_LB[0])*factor); + translation_samp.push_back(translation_LB[1]+(translation_UB[1]-translation_LB[1])*factor); + translation_samp.push_back(translation_LB[2]+(translation_UB[2]-translation_LB[2])*factor); + // Update the translation coordinates in the current state control point object + gipt.trnsl=translation_samp; + + // Again 3 vectors defined, lower bound, upper bound, and sampled, this time for rotation coordinates + vector rotation_LB=cpts[iindex-1].rot; + vector rotation_UB=cpts[iindex].rot; + vector rotation_samp; + + // Now set rotations. Current coordinates computed as before (lowerbound+n% of interval length), but now we also convert degrees to radians + // These must be done case by case (since order of arguments matters) + if (cpts[iindex].rot.size()==2) { + rotation_samp.push_back((rotation_LB[0]+(rotation_UB[0]-rotation_LB[0])*factor)*(M_PI/180)); + rotation_samp.push_back((rotation_LB[1]+(rotation_UB[1]-rotation_LB[1])*factor)*(M_PI/180)); + } + else if (cpts[iindex].rot.size()==3) { + rotation_samp.push_back((rotation_LB[0]+(rotation_UB[0]-rotation_LB[0])*factor)*(M_PI/180)); + rotation_samp.push_back((rotation_LB[1]+(rotation_UB[1]-rotation_LB[1])*factor)*(M_PI/180)); + rotation_samp.push_back((rotation_LB[2]+(rotation_UB[2]-rotation_LB[2])*factor)*(M_PI/180)); + } + else { + egsWarning("EGS_DynamicShape: Invalid number of rotation parameters\n"); + } + gipt.rot=rotation_samp; + + // Update the rotation coordinates in the current state control point object + return 0; + }; + +} // extern "C" + diff --git a/HEN_HOUSE/egs++/shapes/egs_dynamic_shape/egs_dynamic_shape.h b/HEN_HOUSE/egs++/shapes/egs_dynamic_shape/egs_dynamic_shape.h new file mode 100644 index 000000000..0f8daa4a0 --- /dev/null +++ b/HEN_HOUSE/egs++/shapes/egs_dynamic_shape/egs_dynamic_shape.h @@ -0,0 +1,262 @@ +/* +############################################################################### +# +# EGSnrc egs++ dynamic shape 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: Reid Townson +# +# Contributors: +# +############################################################################### +*/ + + +/*! \file egs_dynamic_shape.h + * \brief A dynamic shape + * \RT + */ + +#ifndef EGS_DYNAMIC_SHAPE_ +#define EGS_DYNAMIC_SHAPE_ + +#include "egs_shapes.h" +#include "egs_rndm.h" +#include "egs_application.h" + +#ifdef WIN32 + + #ifdef BUILD_DYNAMIC_SHAPE_DLL + #define EGS_DYNAMIC_SHAPE_EXPORT __declspec(dllexport) + #else + #define EGS_DYNAMIC_SHAPE_EXPORT __declspec(dllimport) + #endif + #define EGS_DYNAMIC_SHAPE_LOCAL + +#else + + #ifdef HAVE_VISIBILITY + #define EGS_DYNAMIC_SHAPE_EXPORT __attribute__ ((visibility ("default"))) + #define EGS_DYNAMIC_SHAPE_LOCAL __attribute__ ((visibility ("hidden"))) + #else + #define EGS_DYNAMIC_SHAPE_EXPORT + #define EGS_DYNAMIC_SHAPE_LOCAL + #endif + +#endif + +/*! \brief An dynamic shape + +\ingroup Shapes + +An dynamic shape is a shape that +takes a random point from another shape and then +applies a transformation, using a time sampling and interpolating between +control points. + +An dynamic shape is defined using +\verbatim +:start shape: + library = egs_dynamic_shape + :start shape: + definition of the shape to be 'dynamic' + :stop shape: + :start motion: # units of cm and degrees + control point = time(1) xtrans(1) ytrans(1) ztrans(1) xrot(1) yrot(1) zrot(1) + control point = time(2) xtrans(2) ytrans(2) ztrans(2) xrot(2) yrot(2) zrot(2) + . + . + . + control point = time(N) xtrans(N) ytrans(N) ztrans(N) xrot(N) yrot(N) zrot(N) + :stop motion: +:stop shape: +\endverbatim + +Control points must be defined such that time(i+1)>=time(i), where time(i) +is the value of time for control point i. The time(i) are automatically +normalized by time(N), where N is the number of control points. + +A translation from the starting position of the shape is applied according to +x, y and z. A rotation follows the same rotation technique as in +EGS_AffineTransform, using the rotation input parameter for 2 or 3 values. +Angles are in degrees and translations in cm. + +Continuous, dynamic motion between control points is simulated by choosing a random +number, R, on (0,1] and, for time(i)0.0, in the case where a user desires to eliminate particles +associated with a range of time values, but there will be a lot of +warning messages. + +*/ +class EGS_DYNAMIC_SHAPE_EXPORT EGS_DynamicShape : public EGS_BaseShape { + +public: + + /*! + * \brief Constructor for EGS_DynamicShape + * \param Shape Base shape to be made dynamic + * \param dyninp Input containing dynamic shape specifications + * \param Name Name of the dynamic shape + * \param f EGS_ObjectFactory pointer + */ + EGS_DynamicShape(EGS_BaseShape *Shape, EGS_Input *dyninp, const string &Name="",EGS_ObjectFactory *f=0) : + EGS_BaseShape(Name,f), shape(Shape) { + if (shape) { + shape->ref(); + otype = "dynamic "; + otype += shape->getObjectType(); + } + else { + otype = "Invalid DynamicShape"; + } + + // Build dynamic shape is where many of the geometry attributes (including control points) are extracted + buildDynamicShape(dyninp); + + if (cpts.size() < 2) { + egsWarning("DynamicShape: not enough or missing control points.\n"); + } + else { + if (cpts[0].time > 0.0) { + egsWarning("DynamicShape: time index of control point 1 > 0.0. This will generate many warning messages.\n"); + } + int npts = cpts.size(); + for (int i=0; i0 && cpts[i].time < cpts[i-1].time-epsilon) { + egsWarning("DynamicShape: time index of control point %i < time index of control point %i\n",i,i-1); + } + if (cpts[i].time<0.0) { + egsWarning("DynamicShape: time index of control point %i < 0.0\n",i); + } + } + + // Normalize time values + for (int i=0; igetPoint(rndm)); + return v; + }; + + /*! \brief Returns a random 3D vector. + * + * Uses the function getPoint() to pick a random position and + * then applies the affine transformation attached to the shape before + * returning it. + */ + EGS_Vector getRandomPoint(EGS_RandomGenerator *rndm) { + getNextShapePosition(rndm); + return shape->getRandomPoint(rndm); + }; + + /*! + * \brief Structure representing a control point for dynamic motion + */ + struct EGS_ControlPoint { + EGS_Float time; //!< Time index for control point + vector trnsl; //!< Vector specifying x, y, z translation + vector rot; //!< Rotation vector + }; + + /*! + * \brief Get the direction of the point source for a given position + * \param Xo Position vector + * \param rndm Random number generator + * \param u Direction vector + * \param wt Weight + */ + void getPointSourceDirection(const EGS_Vector &Xo, + EGS_RandomGenerator *rndm, EGS_Vector &u, EGS_Float &wt) { + if (shape->supportsDirectionMethod()) { + getNextShapePosition(rndm); + shape->getPointSourceDirection(Xo, rndm, u, wt); + } + } + +protected: + + EGS_BaseShape *shape; //!< Base shape made dynamic + + vector cpts; //!< Control points + + int ncpts; //!< Number of control points + + EGS_Float ptime; //!< Time index corresponding to particle + + /*! + * \brief Get the next state of the dynamic shape + * \param rndm Random number generator + */ + void getNextShapePosition(EGS_RandomGenerator *rndm); + + /*! + * \brief Determine whether the simulation geometry contains a dynamic shape + * \param hasdynamic Boolean indicating if the simulation geometry contains a dynamic shape + */ + void containsDynamic(bool &hasdynamic) { + hasdynamic = true; + } + + /*! + * \brief Check if the shape supports the direction method + * \return Boolean indicating if the shape supports the direction method + */ + bool supportsDirectionMethod() const { + return shape->supportsDirectionMethod(); + } + + /*! + * \brief Extract coordinates for the next dynamic shape position + * \param rand Random number for time sampling + * \param gipt EGS_ControlPoint structure to store the coordinates + * \return 0 if successful, otherwise 1 + */ + int getCoord(EGS_Float rand, EGS_ControlPoint &gipt); + + /*! + * \brief Build the dynamic shape using input specifications + * \param dyninp Input containing dynamic shape specifications + */ + void buildDynamicShape(EGS_Input *dyninp); + +}; + +#endif diff --git a/HEN_HOUSE/egs++/sources/egs_beam_source/egs_beam_source.cpp b/HEN_HOUSE/egs++/sources/egs_beam_source/egs_beam_source.cpp index 813489d2e..41dbc44f0 100644 --- a/HEN_HOUSE/egs++/sources/egs_beam_source/egs_beam_source.cpp +++ b/HEN_HOUSE/egs++/sources/egs_beam_source/egs_beam_source.cpp @@ -27,6 +27,7 @@ # Frederic Tessier # Reid Townson # Ernesto Mainegra-Hing +# Alexandre Demelo # ############################################################################### */ @@ -58,7 +59,7 @@ EGS_BeamSource::EGS_BeamSource(EGS_Input *input, EGS_ObjectFactory *f) : i_reuse_photon = 0; i_reuse_electron = 0; is_valid = false; - mu_stored = false; + time_stored = false; lib = 0; Xmin = -veryFar; Xmax = veryFar; @@ -194,15 +195,15 @@ EGS_BeamSource::EGS_BeamSource(EGS_Input *input, EGS_ObjectFactory *f) : n_reuse_electron = ntmp; } - //now see if Mu is returned from the BEAM source - //use motionsample function, tmu will be -1 if not provided by BEAM - //have to save the values read here to use as the first particle in the - //simulation - motionsample(&tei,&txi,&tyi,&tzi,&tui,&tvi,&twi,&twti,&tqi,&tlatchi,&counti,&tiphati,&tmui); - if (!mu_stored) { - if (tmui >= 0.0 && tmui <= 1.0) { - egsInformation("EGS_BeamSource:: Mu index passed from this source.\n"); - mu_stored = true; + // now see if time is returned from the BEAM source use motionsample + // function, ttime will be -1 if not provided by BEAM have to save the + // values read here to use as the first particle in the simulation + motionsample(&tei,&txi,&tyi,&tzi,&tui,&tvi,&twi,&twti,&tqi,&tlatchi,&counti,&tiphati,&ttimei); + + if (!time_stored) { + if (ttimei >= 0.0 && ttimei <= 1.0) { + egsInformation("EGS_BeamSource:: Time index passed from this source.\n"); + time_stored = true; } } @@ -224,7 +225,7 @@ EGS_I64 EGS_BeamSource::getNextParticle(EGS_RandomGenerator *, int &q, wt = wt_save; x = x_save; u = u_save; - mu = mu_save; + time = time_save; ++i_reuse_photon; return count; } @@ -235,11 +236,11 @@ EGS_I64 EGS_BeamSource::getNextParticle(EGS_RandomGenerator *, int &q, wt = wt_save; x = x_save; u = u_save; - mu = mu_save; + time = time_save; ++i_reuse_electron; return count; } - EGS_Float te,tx,ty,tz,tu,tv,tw,twt,tmu; + EGS_Float te,tx,ty,tz,tu,tv,tw,twt,ttime; int tq,tlatch,tiphat; bool ok; do { @@ -257,20 +258,20 @@ EGS_I64 EGS_BeamSource::getNextParticle(EGS_RandomGenerator *, int &q, tlatch=tlatchi; count=counti; tiphat=tiphati; - tmu=tmui; + ttime=ttimei; use_iparticle=false; } else { - motionsample(&te,&tx,&ty,&tz,&tu,&tv,&tw,&twt,&tq,&tlatch,&count,&tiphat,&tmu); + motionsample(&te,&tx,&ty,&tz,&tu,&tv,&tw,&twt,&tq,&tlatch,&count,&tiphat,&ttime); //sample(&te,&tx,&ty,&tz,&tu,&tv,&tw,&twt,&tq,&tlatch,&count,&tiphat); //egsInformation("EGS_BeamSource::getNextParticle: Got E=%g q=%d wt=%g" // " x=(%g,%g,%g) latch=%d count=%lld\n",te,tq,twt,tx,ty,tz, // tlatch,count); - if (mu_stored && (tmu < 0. || tmu > 1.0)) { + if (time_stored && (ttime < 0. || ttime > 1.0)) { //something's wrong - egsWarning("EGS_BeamSource::getNextParticle: Mu index is stored in this source but mu returned < 0\n"); - egsWarning("Will no longer read mu.\n"); - mu_stored = false; + egsWarning("EGS_BeamSource::getNextParticle: Time index is stored in this source but time returned < 0\n"); + egsWarning("Will no longer read time.\n"); + time_stored = false; } } if (tq) { @@ -314,7 +315,7 @@ EGS_I64 EGS_BeamSource::getNextParticle(EGS_RandomGenerator *, int &q, wt = twt; x = EGS_Vector(tx,ty,tz); u = EGS_Vector(tu,tv,tw); - mu = tmu; + time = ttime; if (save_it) { q_save = tq; latch_save = 0; @@ -322,8 +323,19 @@ EGS_I64 EGS_BeamSource::getNextParticle(EGS_RandomGenerator *, int &q, wt_save = twt; x_save = x; u_save = u; - mu_save = mu; + time_save = time; + } + + if (time_stored) { + setTimeIndex(time); + /* this is setting the time index using the base source set call. We get + * rid of the local getTimeIndex function and it should allow for saving + * time in the base source like all the other sources do */ } + else { + setTimeIndex(-1); + } + return count; } diff --git a/HEN_HOUSE/egs++/sources/egs_beam_source/egs_beam_source.h b/HEN_HOUSE/egs++/sources/egs_beam_source/egs_beam_source.h index 79b0a0941..a77c973cb 100644 --- a/HEN_HOUSE/egs++/sources/egs_beam_source/egs_beam_source.h +++ b/HEN_HOUSE/egs++/sources/egs_beam_source/egs_beam_source.h @@ -151,14 +151,6 @@ class EGS_BEAM_SOURCE_EXPORT EGS_BeamSource : public EGS_BaseSource { EGS_Float getFluence() const { return count; }; - EGS_Float getMu() { - if (mu_stored) { - return mu; - } - else { - return -1.0; - } - }; bool storeState(ostream &data) const { return egsStoreI64(data,count); }; @@ -197,11 +189,11 @@ class EGS_BEAM_SOURCE_EXPORT EGS_BeamSource : public EGS_BaseSource { //< to synchronize dynamic source with this bool is_valid; - bool mu_stored; //!< true if mu index stored + bool time_stored; //!< true if time index stored string the_file_name; ifstream the_file; EGS_Float Emax; - EGS_Float mu; + EGS_Float time; EGS_I64 count; // filters @@ -211,7 +203,7 @@ class EGS_BEAM_SOURCE_EXPORT EGS_BeamSource : public EGS_BaseSource { // temporary particle storage int q_save, latch_save; - EGS_Float E_save, wt_save, mu_save; + EGS_Float E_save, wt_save, time_save; EGS_Vector x_save, u_save; // reusing particles @@ -220,8 +212,8 @@ class EGS_BEAM_SOURCE_EXPORT EGS_BeamSource : public EGS_BaseSource { // stored info for first particle read in // need this because we now query the data to see - // if mu index is passed by the source - EGS_Float tei,txi,tyi,tzi,tui,tvi,twi,twti,tmui; + // if time index is passed by the source + EGS_Float tei,txi,tyi,tzi,tui,tvi,twi,twti,ttimei; int tqi,tlatchi,tiphati; EGS_I64 counti; bool use_iparticle; // true if we want to use the above data instead diff --git a/HEN_HOUSE/egs++/sources/egs_dynamic_source/egs_dynamic_source.cpp b/HEN_HOUSE/egs++/sources/egs_dynamic_source/egs_dynamic_source.cpp index 114e44a1f..946c6249e 100644 --- a/HEN_HOUSE/egs++/sources/egs_dynamic_source/egs_dynamic_source.cpp +++ b/HEN_HOUSE/egs++/sources/egs_dynamic_source/egs_dynamic_source.cpp @@ -23,7 +23,7 @@ # # Author: Blake Walters, 2017 # -# Contributors: +# Contributors: Alexandre Demelo # ############################################################################### */ @@ -57,7 +57,7 @@ EGS_DynamicSource::EGS_DynamicSource(EGS_Input *input, } } //now read inputs relevant to dynamic source - //see if user wants to synchronize source with mu read from + //see if user wants to synchronize source with time read from //iaea phsp or beam simulation source vector sync_options; sync_options.push_back("no"); @@ -80,25 +80,39 @@ EGS_DynamicSource::EGS_DynamicSource(EGS_Input *input, int err; int icpts=1; itos << icpts; - string sstring = "control point " + itos.str(); - while (!(err = dyninp->getInput(sstring,point))) { + + string inputTag = "control point"; + string inputTag_backCompat = "control point " + itos.str(); + EGS_Input *currentInput; + + while (true) { + currentInput = dyninp->takeInputItem(inputTag); + if (!currentInput || currentInput->getInput(inputTag, point)) { + currentInput = dyninp->takeInputItem(inputTag_backCompat); + if (!currentInput || currentInput->getInput(inputTag_backCompat, point)) { + delete currentInput; + break; + } + } + delete currentInput; + if (point.size()!=8) { egsWarning("EGS_DynamicSource: control point %i does not specify 8 values.\n",icpts); valid = false; } else { - if (ncpts>0 && point[7] < cpts[ncpts-1].mu) { - egsWarning("EGS_DynamicSource: mu index of control point %i < mu index of control point %i\n",icpts,ncpts); + if (ncpts>0 && point[7] < cpts[ncpts-1].time) { + egsWarning("EGS_DynamicSource: time index of control point %i < time index of control point %i\n",icpts,ncpts); valid = false; } else if (point[7] < 0.) { - egsWarning("EGS_DynamicSource: mu index of control point %i < 0.0\n",icpts); + egsWarning("EGS_DynamicSource: time index of control point %i < 0.0\n",icpts); valid = false; } else { ncpts++; if (ncpts ==1 && point[7] > 0.0) { - egsWarning("EGS_DynamicSource: mu index of control point 1 > 0.0. This will generate many warning messages.\n"); + egsWarning("EGS_DynamicSource: time index of control point 1 > 0.0. This will generate many warning messages.\n"); } //set cpt values cpt.iso = EGS_Vector(point[0],point[1],point[2]); @@ -106,13 +120,13 @@ EGS_DynamicSource::EGS_DynamicSource(EGS_Input *input, cpt.theta = point[4]; cpt.phi = point[5]; cpt.phicol = point[6]; - cpt.mu = point[7]; + cpt.time = point[7]; //add it to the vector of control points cpts.push_back(cpt); icpts++; itos.str(""); itos << icpts; - sstring = "control point " + itos.str(); + inputTag_backCompat = "control point " + itos.str(); } } } @@ -120,14 +134,14 @@ EGS_DynamicSource::EGS_DynamicSource(EGS_Input *input, egsWarning("EGS_DynamicSource: not enough or missing control points.\n"); valid = false; } - if (cpts[ncpts-1].mu == 0.0) { - egsWarning("EGS_DynamicSource: mu index of last control point = 0. Something's wrong.\n"); + if (cpts[ncpts-1].time == 0.0) { + egsWarning("EGS_DynamicSource: time index of last control point = 0. Something's wrong.\n"); valid = false; } else { - //normalize mu index to max. value - for (int i=0; igetSourceDescription(); if (sync) { - description += "\n Source will be synched with mu values read in (if available)."; + description += "\n Source will be synched with time values read in (if available)."; } } } @@ -158,7 +172,7 @@ int EGS_DynamicSource::getCoord(EGS_Float rand, EGS_ControlPoint &ipt) { int iindex=0; int i; for (i=0; i=mu(i), where mu(i) -is the value of mu for control point i. The mu(i) are automatically -normalized by mu(N), where N is the number of control points. +Control points must be defined such that time(i+1)>=time(i), where time(i) +is the value of time for control point i. The time(i) are automatically +normalized by time(N), where N is the number of control points. Continuous, dynamic motion between control points is simulated by choosing a random -number, R, on (0,1] and, for mu(i)0.0, in the case where a user desires to eliminate particles -associated with a range of mu values, but there will be a lot of +only makes sense if time(1)=0.0. However, the source can function +with time(1)>0.0, in the case where a user desires to eliminate particles +associated with a range of time values, but there will be a lot of warning messages. A simple example is shown below. This first defines a monoenergetic @@ -125,7 +126,7 @@ distance, dsource, of 100 cm above the isocentre at (0,0,0). Control points 1 and 2 rotate the source clockwise around the Y-axis (phi=0) through theta=0-360 degrees, while control points 3 and 4 rotate the source clockwise around the Z-axis (phi=90 degrees) through -phi=0-360 degrees. Note that mu(2)-mu(1)=mu(4)-mu(3), so the rotations +phi=0-360 degrees. Note that time(2)-time(1)=time(4)-time(3), so the rotations around Z and Y are carried out for an equal number of incident photons. If the source being made to move dynamically supplies its own monitor unit index (iaea_phsp_source and egs_beam_source only), then the dynamic @@ -153,10 +154,10 @@ to "yes". name = my_source source name = my_parallel_source :start motion: - control point 1 = 0 0 0 100 0 0 0 0 - control point 2 = 0 0 0 100 360 0 0 0.5 - control point 3 = 0 0 0 100 90 0 0 0.5 - control point 4 = 0 0 0 100 90 360 0 1.0 + control point = 0 0 0 100 0 0 0 0 + control point = 0 0 0 100 360 0 0 0.5 + control point = 0 0 0 100 90 0 0 0.5 + control point = 0 0 0 100 90 360 0 1.0 :stop motion: :stop source: @@ -179,7 +180,7 @@ class EGS_DYNAMIC_SOURCE_EXPORT EGS_DynamicSource : EGS_Float theta; //angle of rotation about Y-axis EGS_Float phi; //angle of rotation about Z-axis EGS_Float phicol; //angle of rotation in source plane - EGS_Float mu; //monitor unit index for control point + EGS_Float time; //monitor unit index for control point }; /*! \brief Construct a dynamic source using \a Source as the @@ -195,23 +196,23 @@ class EGS_DYNAMIC_SOURCE_EXPORT EGS_DynamicSource : valid = false; } else { - if (cpts[0].mu > 0.0) { - egsWarning("EGS_DynamicSource: mu index of control point 1 > 0.0. This will generate many warning messages.\n"); + if (cpts[0].time > 0.0) { + egsWarning("EGS_DynamicSource: time index of control point 1 > 0.0. This will generate many warning messages.\n"); } int npts = cpts.size(); for (int i=0; i0 && cpts[i].mu < cpts[i-1].mu) { - egsWarning("EGS_DynamicSource: mu index of control point %i < mu index of control point %i\n",i,i-1); + if (i>0 && cpts[i].time < cpts[i-1].time) { + egsWarning("EGS_DynamicSource: time index of control point %i < time index of control point %i\n",i,i-1); valid = false; } - if (cpts[i].mu<0.0) { - egsWarning("EGS_DynamicSource: mu index of control point %i < 0.0\n",i); + if (cpts[i].time<0.0) { + egsWarning("EGS_DynamicSource: time index of control point %i < 0.0\n",i); valid = false; } } - //normalize mu values + //normalize time values for (int i=0; igetNextParticle(rndm,q,latch,E,wt,x,u); if (sync) { - pmu = source->getMu(); - if (pmu<0) { + ptime = source->getTimeIndex(); + if (ptime<0) { egsWarning("EGS_DynamicSource: You have selected synchronization of dynamic source with %s\n",source->getObjectName().c_str()); - egsWarning("However, this source does not return mu values for each particle. Will turn off synchronization.\n"); + egsWarning("However, this source does not return time values for each particle. Will turn off synchronization.\n"); sync = false; } } if (!sync) { - pmu = rndm->getUniform(); + ptime = rndm->getUniform(); } - err = getCoord(pmu,ipt); + setTimeIndex(ptime); //this is added for the storing of time index in a single location. Technically stored as a basesource attribute not dynamic source + err = getCoord(ptime,ipt); } //translate source in Z @@ -267,9 +269,6 @@ class EGS_DYNAMIC_SOURCE_EXPORT EGS_DynamicSource : EGS_Float getFluence() const { return source->getFluence(); }; - EGS_Float getMu() { - return pmu; - }; bool storeState(ostream &data) const { return source->storeState(data); }; @@ -301,12 +300,12 @@ class EGS_DYNAMIC_SOURCE_EXPORT EGS_DynamicSource : bool valid; //is this a valid source - bool sync; //set to true if source motion synched with mu read from + bool sync; //set to true if source motion synched with time read from //iaea phsp or beam simulation source int getCoord(const EGS_Float rand, EGS_ControlPoint &ipt); - EGS_Float pmu; //monitor unit index corresponding to particle + EGS_Float ptime; //time index corresponding to particle //could just be a random number. void setUp(); diff --git a/HEN_HOUSE/egs++/sources/iaea_phsp_source/iaea_phsp_source.cpp b/HEN_HOUSE/egs++/sources/iaea_phsp_source/iaea_phsp_source.cpp index b691b15e6..5219f15b1 100644 --- a/HEN_HOUSE/egs++/sources/iaea_phsp_source/iaea_phsp_source.cpp +++ b/HEN_HOUSE/egs++/sources/iaea_phsp_source/iaea_phsp_source.cpp @@ -25,6 +25,7 @@ # # Contributors: Reid Townson # Hubert Ho +# Alexandre Demelo # ############################################################################### */ @@ -184,16 +185,16 @@ void IAEA_PhspSource::openFile(const string &phsp_file) { break; } } - //now see if mu index is stored in the file + //now see if time index is stored in the file //assume it is the first variable stored in extra_float //after zlast - i_mu=-1; - mu_stored=false; + i_time=-1; + time_stored=false; if (n_extra_floats > i_zlast+1) { if (extrafloat_types[i_zlast+1] == 0) { - i_mu=i_zlast+1; - mu_stored=true; - egsInformation("IAEA_PhspSource::openFile: Mu index included in data in %s.IAEAphsp\n",phsp_file.c_str()); + i_time=i_zlast+1; + time_stored=true; + egsInformation("IAEA_PhspSource::openFile: Time index included in data in %s.IAEAphsp\n",phsp_file.c_str()); } } @@ -340,8 +341,14 @@ EGS_I64 IAEA_PhspSource::getNextParticle(EGS_RandomGenerator *, int &q, if (mode2) { p.zlast = extrafloattemp[i_zlast]; } - if (mu_stored) { - p.mu = extrafloattemp[i_mu]; + if (time_stored) { + p.time = extrafloattemp[i_time]; + setTimeIndex(p.time); + /* this is setting the time index using the base source set call. We get rid of the local getTimeIndex function and + * it should allow for saving time in the base source like all the other sources do */ + } + else { + setTimeIndex(-1); } if (swap_bytes) { egsSwapBytes(&p.q); @@ -356,7 +363,7 @@ EGS_I64 IAEA_PhspSource::getNextParticle(EGS_RandomGenerator *, int &q, egsSwapBytes(&p.u); egsSwapBytes(&p.v); egsSwapBytes(&p.w); - egsSwapBytes(&p.mu); + egsSwapBytes(&p.time); } //do check here because we need the swapped version of nstat if (nstat<0) { diff --git a/HEN_HOUSE/egs++/sources/iaea_phsp_source/iaea_phsp_source.h b/HEN_HOUSE/egs++/sources/iaea_phsp_source/iaea_phsp_source.h index bc94f0356..2c9d2790c 100644 --- a/HEN_HOUSE/egs++/sources/iaea_phsp_source/iaea_phsp_source.h +++ b/HEN_HOUSE/egs++/sources/iaea_phsp_source/iaea_phsp_source.h @@ -160,14 +160,6 @@ class IAEA_PHSP_SOURCE_EXPORT IAEA_PhspSource : public EGS_BaseSource { double aux = ((double) Nread)/((double) Nparticle); return Pinc*aux; }; - EGS_Float getMu() { - if (mu_stored) { - return p.mu; - } - else { - return -1.0; - } - }; bool storeState(ostream &data) const { data << endl; bool res = egsStoreI64(data,Nread); @@ -251,7 +243,7 @@ class IAEA_PHSP_SOURCE_EXPORT IAEA_PhspSource : public EGS_BaseSource { string the_file_name; //!< The phase-space file name bool mode2; //!< \c true, if a MODE2 file (i.e. storing Zlast) bool latch_stored; //!< true if LATCH is stored in data - bool mu_stored; //!< true if mu index stored + bool time_stored; //!< true if time index stored bool swap_bytes; /*!< \c true, if phase-space file was generated on a CPU with different endianness */ EGS_Float Emax, //!< Maximum k.e. (obtained from the phsp file) @@ -274,7 +266,7 @@ class IAEA_PHSP_SOURCE_EXPORT IAEA_PhspSource : public EGS_BaseSource { int n_extra_floats; //!< no. of extra floats stored in phsp file int n_extra_longs; //!< no. of extra longs stored in phsp file int i_zlast; //!< index of zlast in extra_floats array - int i_mu; //!< index of mu index in extra_floats array + int i_time; //!< index of time index in extra_floats array int i_latch; //!< index of latch in extra_floats array bool first; @@ -292,7 +284,7 @@ class IAEA_PHSP_SOURCE_EXPORT IAEA_PhspSource : public EGS_BaseSource { #ifndef SKIP_DOXYGEN struct EGS_LOCAL BeamParticle { int latch, q; - float E, u, v, w, x, y, z, wt, zlast, mu; + float E, u, v, w, x, y, z, wt, zlast, time; }; BeamParticle p; #endif diff --git a/HEN_HOUSE/egs++/view/egs_track_view.cpp b/HEN_HOUSE/egs++/view/egs_track_view.cpp index 02660c6f7..b1fd2289b 100644 --- a/HEN_HOUSE/egs++/view/egs_track_view.cpp +++ b/HEN_HOUSE/egs++/view/egs_track_view.cpp @@ -26,6 +26,7 @@ # Contributors: Iwan Kawrakow # Frederic Tessier # Ernesto Mainegra-Hing +# Alexandre Demelo # ############################################################################### */ @@ -67,7 +68,16 @@ static T *shrink(T *original, size_t old_len, size_t new_len) { const int zero = 0; -EGS_TrackView::EGS_TrackView(const char *filename, vector &ntracks) { +EGS_TrackView::EGS_TrackView(const char *filename, vector &ntracks,vector &timelist_p, vector &timelist_e, vector &timelist_po) { + // timelist vector references added as arguments in order to incorporate + // particle track time slider in egsview. + + /* time lists are references so all changes here change the time list + * vectors in renderworker::loadtracks. The lists then pass from + * renderworker to viewcontrol through signal emissions. tracksloaded from + * renderworker to imagewindow tracksloaded from image window to viewcontrol + */ + // typedefs to keep things short typedef EGS_ParticleTrack::Vertex Vert; typedef EGS_ParticleTrack::ParticleInfo PInfo; @@ -93,6 +103,7 @@ EGS_TrackView::EGS_TrackView(const char *filename, vector &ntracks) { data.seekg(0, ios::beg); data.read((char *)&tot_tracks, sizeof(int)); + // very conservative sanity check to avoid a huge allocation if (tot_tracks * sizeof(Vert) > size - sizeof(int)) { egsInformation("%s: No tracks loaded: %d tracks can't fit in %d-byte file '%s'\n", @@ -102,6 +113,13 @@ EGS_TrackView::EGS_TrackView(const char *filename, vector &ntracks) { } egsInformation("%s: Reading %d tracks from '%s' ...\n", func_name, tot_tracks, filename); + string readfile = string(filename); + string extension; + size_t k = readfile.rfind('.', readfile.length()); + if (k != string::npos) { + extension = readfile.substr(k+1, readfile.length() - k); + } + tmp_buffer = new char[size-sizeof(int)]; // May want to look into memory mapping, but only if access patterns/OS sets matter @@ -129,12 +147,29 @@ EGS_TrackView::EGS_TrackView(const char *filename, vector &ntracks) { char *loc = tmp_buffer; for (int i=0; i &ntracks) { count_vert[type] += nverts; tmp_index[i] = loc; - int skip = sizeof(int) + sizeof(PInfo); + loc += skip; for (int k=0; ke; @@ -175,11 +210,22 @@ EGS_TrackView::EGS_TrackView(const char *filename, vector &ntracks) { // Compression routine! for (int i=0; i &ntracks) { m_index[type][i_off] = m_off; ind_rcnt[type]++; mem_rcnt[type] += nverts; + if (extension=="syncptracks") { + // update timelist vectors such that the time indices of the + // compressed particle list is available from the viewcontrol + // class lists updated here to ensure the time index vector + // indices match those of the compressed particle list. + if (type==0) { + timelist_p.push_back(timeindex); + } + else if (type==1) { + timelist_e.push_back(timeindex); + } + else if (type==2) { + timelist_po.push_back(timeindex); + } + + } } } // Sentinel at the end; get lengths @@ -240,8 +302,11 @@ EGS_TrackView::EGS_TrackView(const char *filename, vector &ntracks) { ind_rcnt[0]+ind_rcnt[1]+ind_rcnt[2], ind_rcnt[0], ind_rcnt[1], ind_rcnt[2]); m_failed = false; + data.close(); + } + EGS_TrackView::~EGS_TrackView() {} bool EGS_TrackView::renderTracks(int nx, int ny, EGS_Vector *image, diff --git a/HEN_HOUSE/egs++/view/egs_track_view.h b/HEN_HOUSE/egs++/view/egs_track_view.h index 6cd0e3565..a537e7770 100644 --- a/HEN_HOUSE/egs++/view/egs_track_view.h +++ b/HEN_HOUSE/egs++/view/egs_track_view.h @@ -25,6 +25,7 @@ # # Contributors: Iwan Kawrakow # Frederic Tessier +# Alexandre Demelo # ############################################################################### */ @@ -38,6 +39,7 @@ #include "egs_transformations.h" #include "stddef.h" #include "egs_particle_track.h" +#include class EGS_Matrix : private EGS_RotationMatrix { public: @@ -92,11 +94,11 @@ class EGS_Matrix : private EGS_RotationMatrix { class EGS_TrackView { public: - - EGS_TrackView(const char *filename, vector &ntracks); + EGS_TrackView(const char *filename, vector &ntracks, vector &timelist_p, vector &timelist_e, vector &timelist_po); ~EGS_TrackView(); + bool renderTracks(int nx, int ny, EGS_Vector *image, EGS_ClippingPlane **planes, const int n_planes, int *abort_location=NULL); @@ -123,7 +125,9 @@ class EGS_TrackView { trackIndices = trackInd; } + protected: + void renderTrack(EGS_ParticleTrack::Vertex *const vs, int len, EGS_Float color, int nx, int ny, EGS_Vector *image); // High-level camera description @@ -146,6 +150,7 @@ class EGS_TrackView { vector trackIndices; + EGS_ClippingPlane m_planes[14]; // Clipping planes. 0-3 are for the viewport int nplanes; // number of planes used diff --git a/HEN_HOUSE/egs++/view/egs_visualizer.cpp b/HEN_HOUSE/egs++/view/egs_visualizer.cpp index 662a269e8..3e89bb49c 100644 --- a/HEN_HOUSE/egs++/view/egs_visualizer.cpp +++ b/HEN_HOUSE/egs++/view/egs_visualizer.cpp @@ -24,6 +24,7 @@ # Author: Iwan Kawrakow, 2005 # # Contributors: Frederic Tessier +# Alexandre Demelo # ############################################################################### */ @@ -42,12 +43,12 @@ class EGS_PrivateVisualizer { nlight(0), ntot(0), nmat(0), nclip(0), nclip_t(0), m_tracks(NULL) {}; ~EGS_PrivateVisualizer(); - vector loadTracksData(const char *fname) { + vector loadTracksData(const char *fname, vector &timelist_p, vector &timelist_e, vector &timelist_po) { if (m_tracks) { delete m_tracks; } vector ntracks; - m_tracks = new EGS_TrackView(fname, ntracks); + m_tracks = new EGS_TrackView(fname, ntracks, timelist_p, timelist_e, timelist_po); return ntracks; } void setProjection(const EGS_Vector &camera_pos, @@ -134,6 +135,7 @@ class EGS_PrivateVisualizer { // get the color for the screen point x. EGS_Vector getColor(const EGS_Vector &x, EGS_BaseGeometry *g, const EGS_Float track_distance, const EGS_Float track_alpha, bool debug=false, EGS_Vector bCol=EGS_Vector(1,1,1)); + void getRegions(const EGS_Vector &x, EGS_BaseGeometry *g, int *regions, EGS_Vector *colors, int maxreg, EGS_Vector &hitCoord, const unordered_map &score, EGS_Float &hitScore); void getFirstHit(const EGS_Vector &x, EGS_BaseGeometry *g, EGS_Vector &hitCoord); @@ -206,10 +208,10 @@ EGS_GeometryVisualizer::~EGS_GeometryVisualizer() { delete p; } -vector EGS_GeometryVisualizer::loadTracksData(const char *fname) { +vector EGS_GeometryVisualizer::loadTracksData(const char *fname, vector &timelist_p, vector &timelist_e, vector &timelist_po) { vector ntracks; if (p) { - ntracks = p->loadTracksData(fname); + ntracks = p->loadTracksData(fname,timelist_p, timelist_e, timelist_po); } return ntracks; } diff --git a/HEN_HOUSE/egs++/view/egs_visualizer.h b/HEN_HOUSE/egs++/view/egs_visualizer.h index f7f257f93..659db6bc6 100644 --- a/HEN_HOUSE/egs++/view/egs_visualizer.h +++ b/HEN_HOUSE/egs++/view/egs_visualizer.h @@ -24,6 +24,7 @@ # Author: Iwan Kawrakow, 2005 # # Contributors: Frederic Tessier +# Alexandre Demelo # ############################################################################### */ @@ -88,7 +89,7 @@ class EGS_GeometryVisualizer { EGS_GeometryVisualizer(); ~EGS_GeometryVisualizer(); - vector loadTracksData(const char *fname); + vector loadTracksData(const char *fname, vector &timelist_p, vector &timelist_e, vector &timelist_po); void setProjection(const EGS_Vector &camera_pos, const EGS_Vector &camera_look_at, EGS_Float distance, @@ -130,7 +131,6 @@ class EGS_GeometryVisualizer { // region picking void regionPick(int x, int y); - #ifdef HAVE_PNG bool makePngImage(EGS_BaseGeometry *, int xsize, int ysize, const char *fname); diff --git a/HEN_HOUSE/egs++/view/image_window.cpp b/HEN_HOUSE/egs++/view/image_window.cpp index 8c6ed5fe1..dae6922ac 100644 --- a/HEN_HOUSE/egs++/view/image_window.cpp +++ b/HEN_HOUSE/egs++/view/image_window.cpp @@ -26,6 +26,7 @@ # Contributors: Frederic Tessier # Manuel Stoeckl # Reid Townson +# Alexandre Demelo # ############################################################################### */ @@ -87,6 +88,7 @@ ImageWindow::ImageWindow(QWidget *parent, const char *name) : qRegisterMetaType("RenderParameters"); qRegisterMetaType("RenderResults"); qRegisterMetaType>("vector"); + qRegisterMetaType>("vector"); // Initialize render worker and put it in a thread restartWorker(); @@ -225,6 +227,7 @@ void ImageWindow::saveView(EGS_BaseGeometry *geo, int nx, int ny, QString name, void ImageWindow::loadTracks(QString name) { emit requestLoadTracks(name); + } void ImageWindow::stopWorker() { @@ -250,7 +253,7 @@ void ImageWindow::restartWorker() { worker, SLOT(render(EGS_BaseGeometry *,RenderParameters))); connect(this, SIGNAL(requestLoadTracks(QString)), worker, SLOT(loadTracks(QString))); connect(worker, SIGNAL(rendered(RenderResults,RenderParameters)), this, SLOT(drawResults(RenderResults,RenderParameters))); - connect(worker, SIGNAL(tracksLoaded(vector)), this, SLOT(trackResults(vector))); + connect(worker, SIGNAL(tracksLoaded(vector, vector, vector, vector)), this, SLOT(trackResults(vector, vector, vector, vector))); connect(worker, SIGNAL(aborted()), this, SLOT(handleAbort())); thread->start(); renderState = WorkerIdle; @@ -647,8 +650,8 @@ void ImageWindow::drawResults(RenderResults r, RenderParameters q) { repaint(); } -void ImageWindow::trackResults(vector ntracks) { - emit tracksLoaded(ntracks); +void ImageWindow::trackResults(vector ntracks, vector timelist_p, vector timelist_e, vector timelist_po) { + emit tracksLoaded(ntracks, timelist_p, timelist_e, timelist_po); } void ImageWindow::handleAbort() { diff --git a/HEN_HOUSE/egs++/view/image_window.h b/HEN_HOUSE/egs++/view/image_window.h index 38d082fed..837f2d0dd 100644 --- a/HEN_HOUSE/egs++/view/image_window.h +++ b/HEN_HOUSE/egs++/view/image_window.h @@ -25,6 +25,7 @@ # # Contributors: Frederic Tessier # Manuel Stoeckl +# Alexandre Demelo # ############################################################################### */ @@ -90,7 +91,7 @@ public slots: protected slots: void drawResults(RenderResults,RenderParameters); - void trackResults(vector); + void trackResults(vector, vector, vector, vector); void handleAbort(); signals: @@ -106,7 +107,7 @@ protected slots: void leftMouseClick(int x, int y); void leftDoubleClick(EGS_Vector hitCoord); void saveComplete(); - void tracksLoaded(vector); + void tracksLoaded(vector, vector, vector, vector); // for render thread void requestRender(EGS_BaseGeometry *,RenderParameters); diff --git a/HEN_HOUSE/egs++/view/main.cpp b/HEN_HOUSE/egs++/view/main.cpp index 3263198a0..950d05ce8 100644 --- a/HEN_HOUSE/egs++/view/main.cpp +++ b/HEN_HOUSE/egs++/view/main.cpp @@ -24,6 +24,7 @@ # Author: Iwan Kawrakow, 2005 # # Contributors: Frederic Tessier +# Alexandre Demelo # ############################################################################### */ @@ -100,6 +101,7 @@ int main(int argc, char **argv) { QString tracks_file = QString(""); QString config_file = QString(""); + QString extension=QString(""); if (argc >= 3) { QString argv2 = argv[2]; if (argv2.endsWith("ptracks")) { @@ -119,7 +121,17 @@ int main(int argc, char **argv) { } } + // get the extension of the tracks file (either ptracks or syncptracks) if + // tracks file provided in terminal to pass to the viewcontrol + if (tracks_file.endsWith("syncptracks")) { + extension=QString("syncptracks"); + } + else { + extension=QString("ptracks"); + } + w.setTracksFilename(tracks_file); + w.setTracksExtension(extension); // set the tracks extension in the viewcontrol if (!w.loadInput(false)) { return 1; } diff --git a/HEN_HOUSE/egs++/view/renderworker.cpp b/HEN_HOUSE/egs++/view/renderworker.cpp index eda5eb17c..1db48bdc2 100644 --- a/HEN_HOUSE/egs++/view/renderworker.cpp +++ b/HEN_HOUSE/egs++/view/renderworker.cpp @@ -23,6 +23,8 @@ # # Author: Manuel Stoeckl, 2015 # +# Contributors: Alexandre Demelo +# ############################################################################### */ @@ -52,8 +54,12 @@ RenderWorker::~RenderWorker() { } void RenderWorker::loadTracks(QString fileName) { - vector ntracks = vis->loadTracksData(fileName.toUtf8().constData()); - emit tracksLoaded(ntracks); + vector timelist_p, timelist_e, timelist_po; + vector ntracks = vis->loadTracksData(fileName.toUtf8().constData(), timelist_p, timelist_e, timelist_po); + // time lists passed as arguments above are references in all the subsequent + // calls (Geometryvisualizer, then privateVisualizer, then Trackview + // creation which updates time list reference) + emit tracksLoaded(ntracks, timelist_p, timelist_e, timelist_po); } void RenderWorker::drawAxes(const RenderParameters &p) { @@ -340,4 +346,3 @@ struct RenderResults RenderWorker::renderSync(EGS_BaseGeometry *g, struct Render return r; } - diff --git a/HEN_HOUSE/egs++/view/renderworker.h b/HEN_HOUSE/egs++/view/renderworker.h index b65ef573d..6a0a00472 100644 --- a/HEN_HOUSE/egs++/view/renderworker.h +++ b/HEN_HOUSE/egs++/view/renderworker.h @@ -23,6 +23,7 @@ # # Author: Manuel Stoeckl, 2015 # +# Contributors: Alexandre Demelo ############################################################################### */ @@ -135,7 +136,7 @@ public slots: void aborted(); void rendered(struct RenderResults, struct RenderParameters params); - void tracksLoaded(vector ntracks); + void tracksLoaded(vector ntracks, vector timelist_p, vector timelist_e, vector timelist_po); private: void drawAxes(const struct RenderParameters &); diff --git a/HEN_HOUSE/egs++/view/viewcontrol.cpp b/HEN_HOUSE/egs++/view/viewcontrol.cpp index 5d3bd48ab..c5c643b59 100644 --- a/HEN_HOUSE/egs++/view/viewcontrol.cpp +++ b/HEN_HOUSE/egs++/view/viewcontrol.cpp @@ -27,6 +27,7 @@ # Ernesto Mainegra-Hing # Manuel Stoeckl # Reid Townson +# Alexandre Demelo # ############################################################################### */ @@ -35,6 +36,7 @@ #include "saveimage.h" #include "clippingplanes.h" #include "viewcontrol.h" +#include "egs_track_view.h" #include "egs_libconfig.h" #include "egs_functions.h" @@ -44,6 +46,7 @@ #include "egs_input.h" #include "egs_ausgab_object.h" #include "ausgab_objects/egs_dose_scoring/egs_dose_scoring.h" +#include "egs_particle_track.h" #include #include @@ -62,6 +65,10 @@ #include #include #include +#include +#include +#include + using namespace std; #ifndef M_PI @@ -125,7 +132,6 @@ GeometryViewControl::GeometryViewControl(QWidget *parent, const char *name) if (showPhotonTracks || showElectronTracks || showPositronTracks) { showTracks = true; } - // camera orientation vectors (same as the screen vectors) camera_v1 = screen_v1; camera_v2 = screen_v2; @@ -164,7 +170,10 @@ GeometryViewControl::GeometryViewControl(QWidget *parent, const char *name) connect(gview, SIGNAL(putCameraOnAxis(char)), this, SLOT(cameraOnAxis(char))); connect(gview, SIGNAL(leftMouseClick(int,int)), this, SLOT(reportViewSettings(int,int))); connect(gview, SIGNAL(leftDoubleClick(EGS_Vector)), this, SLOT(setRotationPoint(EGS_Vector))); - connect(gview, SIGNAL(tracksLoaded(vector)), this, SLOT(updateTracks(vector))); + + // EGS_Float vectors added to signal & slot to allow the track loading + // method in trackview to set the time index lists in this class + connect(gview, SIGNAL(tracksLoaded(vector, vector, vector, vector)), this, SLOT(updateTracks(vector, vector, vector, vector))); save_image = new SaveImage(this,"save image"); @@ -177,6 +186,9 @@ GeometryViewControl::GeometryViewControl(QWidget *parent, const char *name) // set the widget to show near the left-upper corner of the screen move(QPoint(25,25)); + + //set the play button active boolean to false + isPlaying=false; } GeometryViewControl::~GeometryViewControl() { @@ -221,7 +233,6 @@ bool GeometryViewControl::loadInput(bool reloading, EGS_BaseGeometry *simGeom) { #ifdef VIEW_DEBUG egsWarning("In loadInput(), reloading is %d\n",reloading); #endif - // check that the file (still) exists QFile file(filename); if (!file.exists()) { @@ -230,15 +241,12 @@ bool GeometryViewControl::loadInput(bool reloading, EGS_BaseGeometry *simGeom) { } QFileInfo fileInfo = QFileInfo(file); - // Set the title of the viewer this->setProperty("windowTitle", "View Controls ("+fileInfo.baseName()+")"); gview->setProperty("windowTitle", "egs_view ("+fileInfo.baseName()+")"); - // Clear the current geometry gview->stopWorker(); qApp->processEvents(); - // Delete any previous geometry if (!simGeom) { #ifdef VIEW_DEBUG @@ -250,6 +258,10 @@ bool GeometryViewControl::loadInput(bool reloading, EGS_BaseGeometry *simGeom) { } EGS_BaseGeometry::clearGeometries(); + // solve loading new input file from file tab in egsview. Previously + // name array is never cleared after loading new input file + geometryNames.clear(); + // Delete any previous ausgab objects #ifdef VIEW_DEBUG egsInformation("GeometryViewControl::loadInput: Clearing previous ausgab objects...\n"); @@ -261,14 +273,12 @@ bool GeometryViewControl::loadInput(bool reloading, EGS_BaseGeometry *simGeom) { } } } - // Read the input file #ifdef VIEW_DEBUG egsInformation("GeometryViewControl::loadInput: Reading input file...\n"); #endif EGS_Input input; input.setContentFromFile(filename.toUtf8().constData()); - // Load the new geometry #ifdef VIEW_DEBUG egsInformation("GeometryViewControl::loadInput: Creating the geometry...\n"); @@ -305,7 +315,6 @@ bool GeometryViewControl::loadInput(bool reloading, EGS_BaseGeometry *simGeom) { else { newGeom = simGeom; } - // restart from scratch (copied from main.cpp) EGS_Float xmin = -50, xmax = 50; EGS_Float ymin = -50, ymax = 50; @@ -357,7 +366,6 @@ bool GeometryViewControl::loadInput(bool reloading, EGS_BaseGeometry *simGeom) { } delete vc; } - // Only allow region selection for up to 1k regions int nreg = newGeom->regions(); if (nreg < 1001) { @@ -370,14 +378,11 @@ bool GeometryViewControl::loadInput(bool reloading, EGS_BaseGeometry *simGeom) { egsInformation("Region selection tab has been disabled due to >1000 regions (for performance reasons)\n"); } tabWidget->setTabEnabled(2,allowRegionSelection); - // Get the rendering parameters RenderParameters &rp = gview->pars; rp.allowRegionSelection = allowRegionSelection; rp.trackIndices.assign(6,1); - gview->restartWorker(); - if (!simGeom) { setGeometry(newGeom,user_colors,xmin,xmax,ymin,ymax,zmin,zmax,reloading); origSimGeom = g; @@ -391,14 +396,12 @@ bool GeometryViewControl::loadInput(bool reloading, EGS_BaseGeometry *simGeom) { if (allowRegionSelection) { updateRegionTable(); } - // Set the simulation geometry combobox comboBox_simGeom->clear(); for (unsigned int i=0; iaddItem((geometryNames[i] + " (" + g->getGeometry(geometryNames[i])->getType() + ")").c_str(), geometryNames[i].c_str()); } comboBox_simGeom->setCurrentIndex(comboBox_simGeom->findData(g->getName().c_str())); - // Load ausgab objects from the input file if (!simGeom) { #ifdef VIEW_DEBUG @@ -1082,6 +1085,10 @@ void GeometryViewControl::setTracksFilename(QString str) { filename_tracks = str; } +void GeometryViewControl::setTracksExtension(QString str) { + tracks_extension = str; +} + void GeometryViewControl::loadDose() { // Prompt the user to select a previous config file QFileInfo inputFileInfo = QFileInfo(filename); @@ -1890,20 +1897,44 @@ void GeometryViewControl::loadTracksDialog() { egsWarning("In loadTracksDialog()\n"); #endif QFileInfo inputFileInfo = QFileInfo(filename); - filename_tracks = QFileDialog::getOpenFileName(this, "Select particle tracks file", inputFileInfo.canonicalPath(), "*.ptracks"); + filename_tracks = QFileDialog::getOpenFileName(this, "Select particle tracks file", inputFileInfo.canonicalPath(), "*ptracks"); + // tracks_extension is set based on selected tracks file. syncptracks + // contain track time indices, ptracks do not + if (filename_tracks.endsWith("syncptracks")) { + tracks_extension=QString("syncptracks"); + if (!hasDynamic) { + // if hasdynamic is not yet true (no dynamic geometry) and extension + // is syncptracks then check the time indices in the file. check + // done through hasValidTime function + hasDynamic=hasValidTime(); + } + } + else { + tracks_extension=QString("ptracks"); + } if (filename_tracks.isEmpty()) { return; } + // run timeObjectVisibility to either make visible or hide time index + // related objects depending on input file and tracks file + timeObjectVisibility(); gview->loadTracks(filename_tracks); } -void GeometryViewControl::updateTracks(vector ntracks) { +void GeometryViewControl::updateTracks(vector ntracks, vector timeindexlist_p, vector timeindexlist_e, vector timeindexlist_po) { if (ntracks.size() != 3) { return; } + // vectors containing the sorted list of time indices corresponding to the + // compressed particle tracks list is saved + timelist_p=timeindexlist_p; + timelist_e=timeindexlist_e; + timelist_po=timeindexlist_po; + + #ifdef VIEW_DEBUG egsWarning("In updateTracks(%d %d %d)\n",ntracks[0], ntracks[1], ntracks[2]); #endif @@ -1926,7 +1957,9 @@ void GeometryViewControl::updateTracks(vector ntracks) { spin_tmine->setValue(1); spin_tminpo->setValue(1); - updateView(); + // Update the time window value + // This includes an updateView() call + slideTime(); } void GeometryViewControl::viewAllMaterials() { @@ -2008,9 +2041,27 @@ int GeometryViewControl::setGeometry( } g = geom; + // check the geometry and the tracks file to determine whether the time + // index objects should be made visible (i.e., setting hasdynamic) + hasDynamic=false; + + // loop through the different layers of the geometry and makes hasdynamic + // true if a dynamic geometry is found. This is independent of tracks file + // type, such that visualizing geometry motion is possible even with a + // ptracks file + g->containsDynamic(hasDynamic); if (!filename_tracks.isEmpty()) { gview->loadTracks(filename_tracks); + // if hasdynamic is not yet true (no dynamic geometry) and extension is + // syncptracks then check the time indices in the file. check done + // through hasValidTime function + if (!hasDynamic && tracks_extension=="syncptracks") { + hasDynamic = hasValidTime(); + } } + // run timeObjectVisibility to either make visible or hide time index + // related objects depending on input file and tracks file + timeObjectVisibility(); #ifdef VIEW_DEBUG qDebug("Got %d user defined colors",int(ucolors.size())); @@ -2394,8 +2445,8 @@ int GeometryViewControl::setGeometry( } void GeometryViewControl::updateView(bool transform) { - // transfer RenderParameters &rp = gview->pars; + rp.axesmax = axesmax; rp.camera = camera; rp.camera_v1 = camera_v1; @@ -2716,6 +2767,284 @@ void GeometryViewControl::endTransformation() { gview->endTransformation(); } +// Time index visual elements methods +void GeometryViewControl::playTime() { + if (isPlaying) { + button_timeplay->setText("Play"); + isPlaying=false; + } + else { + button_timeplay->setText("Pause"); + isPlaying=true; + } + int sliderpos=slider_timeindex->sliderPosition(); + + if (sliderpos==spin_numTimeSteps->value()-1) { + sliderpos=0; + } + // this function controls the play button, and allows for the simulation to + // be automatically played out sequentially in time. + for (int i= sliderpos; ivalue();) { + // the simulation plays through 1000 discrete time points (equivalent to + // possible slider steps) from 0.000 to 0.999 in 0.0001 increments + if (!isPlaying) { + break; + } + EGS_Float currtime = i/(float)spin_numTimeSteps->value(); + + // update time index display/input box. The signals are blocked as it + // would lead to an infinite loop between the slider and the time index + // spin box + spin_timeindex->blockSignals(true); + spin_timeindex->setValue(currtime); + spin_timeindex->blockSignals(false); + + // the time index slider's value is updated (as is its position since + // tracking is enabled). Changing the slider value emits a + // 'valuechanged' signal which automatically invokes the slidertime + // function + slider_timeindex->setValue(i); + + // update the Qt display between discrete steps + qApp->processEvents(); + + i+=1; + + // wait 10 milliseconds before making the next step, otherwise the + // motion would be difficult to follow + std::this_thread::sleep_for(std::chrono::milliseconds(10)); + } + button_timeplay->setText("Play"); + isPlaying=false; +} + +void GeometryViewControl::resetTime() { + /* this function controls the reset button, and allows for the time index + * settings to be returned to their initial state. need only reset the slider + * and time window values to their original state. These will automatically + * call the particleslider() and slidetime() methods respectively due to + * 'valuechanged' signal emissions. All other elements (time index spinbox + * and the particle bounds) are returned to their initial states through + * these */ + slider_timeindex->setValue(0); + + // Default the time window to 1% of the simulation time + spin_timewindow->setValue(0.01); +} + +void GeometryViewControl::spinTime() { + // this function controls the time index spinbox, which is used both as a + // way to display the current index, and to input and go to a particular + // index first, get the new time index, and from it calculate the equivalent + // position on the slider (multiply by 1000 to get integer between 0 and + // 999) + EGS_Float slidertime=spin_timeindex->value(); + int sliderpos=(int)(slidertime*spin_numTimeSteps->value()); + + // update the time index slider. The signals are blocked as it would lead + // to an infinite loop between the slider and the time index spin box + slider_timeindex->blockSignals(true); + slider_timeindex->setValue(sliderpos); + slider_timeindex->blockSignals(false); + + // updatePosition is called on our geometry using the new time index, this + // simply sets the geometry transformation as specified by the control + // points + g->updatePosition(slidertime); + + // finally, the particleslider function is called using the new time index + // to update the particle track indices + particleSlider(slidertime); + updateView(); +} + +void GeometryViewControl::slideTime() { + /* this function controls the time index slider and allows for the user to + * control the geometry motion and visible particle tracks by dragging the + * slider to various time indices. Note that this function is connected to + * the slider via the valuechanged signal, such that it will update + * continuously as the slider is dragged, not only at discrete points where + * the slider is released. Note this is also connected to the time window + * spin box. It requires essentially the same code, as it ultimately it only + * needs to call particleslider with the current time index to change the + * particle track indices */ + + // first the slider position is grabbed. The slider ranges from 0 to 999 as + // it cannot have decimal increments + int sliderpos=slider_timeindex->sliderPosition(); + + // the time index corresponding to the position integer is determined by + // dividing the position by 1000 + EGS_Float slidertime = sliderpos/(float)1000; + + // update the time index spinbox. The signals are blocked as it would lead + // to an infinite loop between the slider and the time index spin box + spin_timeindex->blockSignals(true); + spin_timeindex->setValue(slidertime); + spin_timeindex->blockSignals(false); + + // updatePosition is called on our geometry using the new time index, this + // simply sets the geometry transformation as specified by the control + // points + g->updatePosition(slidertime); + + // finally, the particleslider function is called using the new time index + // to update the particle track indices + particleSlider(slidertime); + updateView(); +} + + +void GeometryViewControl::particleSlider(EGS_Float slidertime) { + /* this function is called by the time index slider, the time window box, + * and the time index box, as the final step in updating the display. It is + * responsible for determining which particle tracks are within the range to + * be displayed at a given time index. It considers a range of +/- half the + * provided time window on either side of the slider position and checks the + * sorted list of time indices corresponding to the compressed particle + * tracks list, and determines the start and end index of each particle type + * boxes to be imposed on their min and max spin boxes */ + + /* first, all of these operations can only be performed given a syncptracks + * file is being used as we need time indices to compare too (note a + * syncptracks file is not a given, as one could have access to the time + * index view elements for a dynamic geometry where the time is not + * recorded, however given a syncptracks file we know the time indices are + * valid)*/ + if (tracks_extension=="syncptracks") { + + // obtain size of the time window + EGS_Float t_window = spin_timewindow->value(); + + // a boolean variable hasstart is defined as initially set to false. + // This will track whether a starting index has been assigned + bool hasstart=false; + + // start and end index integers defined + int startindex=1,endindex=1; + + // loop through the photons time index list one at a time + for (int j=0; j(slidertime-(t_window/2))) { + + // if it is, we check whether a starting index has been assigned + if (!hasstart) { + // if the starting index has not been assigned, assign both + // the start and the end index to the current loop index, + // and make hasstart true; the end index is also set here as + // it is possible only one particle is in the range, and the + // end cannot be smaller than the start + startindex=j+1; + endindex=j+2; + hasstart=true; + // since hasstart is now true, the starting index will not + // be updated at any later loop iteration + } + else { + // if the starting index has been assigned, only set the end + // index + endindex=j+2; + } + } + } + // once the loop is over the max and min spin boxes are set (both are 1 + // if no particle in range) (note max must be set first as it sets the + // maximum value of the min spinbox) + spin_tmaxp->blockSignals(true); + spin_tmaxp->setValue(endindex); + spin_tmaxp->blockSignals(false); + spin_tminp->blockSignals(true); + spin_tminp->setValue(startindex); + spin_tminp->blockSignals(false); + + // hasstart bool and indices are reset to initial states + hasstart=false; + startindex=1; + endindex=1; + + // loop through the electrons time index list one at a time + for (int j=0; j(slidertime-(t_window/2))) { + // if it is, we check whether a starting index has been assigned + if (!hasstart) { + // if the starting index has not been assigned, assign both + // the start and the end index to the current loop index, + // and make hasstart true; the end index is also set here as + // it is possible only one particle is in the range, and the + // end cannot be smaller than the start + startindex=j+1; + endindex=j+2; + hasstart=true; + // since hasstart is now true, the starting index will not + // be updated at any later loop iteration + } + else { + // if the starting index has been assigned, only set the end + // index + endindex=j+1; + } + } + } + // once the loop is over the max and min spin boxes are set (both are 1 + // if no particle in range) (note max must be set first as it sets the + // maximum value of the min spinbox) + spin_tmaxe->blockSignals(true); + spin_tmaxe->setValue(endindex); + spin_tmaxe->blockSignals(false); + spin_tmine->blockSignals(true); + spin_tmine->setValue(startindex); + spin_tmine->blockSignals(false); + + // hasstart bool and indices are reset to initial states + hasstart=false; + startindex=1; + endindex=1; + // loop through the positron time index list one at a time + for (int j=0; j(slidertime-(t_window/2))) { + // if it is, we check whether a starting index has been assigned + if (!hasstart) { + // if the starting index has not been assigned, assign both + // the start and the end index to the current loop index, + // and make hasstart true; the end index is also set here as + // it is possible only one particle is in the range, and the + // end cannot be smaller than the start + startindex=j+1; + endindex=j+2; + hasstart=true; + // since hasstart is now true, the starting index will not + // be updated at any later loop iteration + } + else { + // if the starting index has been assigned, only set the end + // index + endindex=j+1; + } + } + } + // once the loop is over the max and min spin boxes are set (both are 1 + // if no particle in range) (note max must be set first as it sets the + // maximum value of the min spinbox) + spin_tmaxpo->blockSignals(true); + spin_tmaxpo->setValue(endindex); + spin_tmaxpo->blockSignals(false); + spin_tminpo->blockSignals(true); + spin_tminpo->setValue(startindex); + spin_tminpo->blockSignals(false); + } + // if the file is not a syncptracks nothing at all occurs in this function + // and the particle tracks are unchanged +} +// end of time index methods// + void GeometryViewControl::updateRegionTable() { if (!g) { return; @@ -2962,3 +3291,71 @@ void GeometryViewControl::setFontSize(int size) { controlsText->setTextCursor(cursor); } +bool GeometryViewControl::hasValidTime() { + /* This function is used to determine whether the time index elements should be shown in the case that no dynamic geometry is present. + * if the tracks file is a .syncptracks file (a condition for callling this function), it will check the first of the time indices written to it. If it is not -1 + * then the time indices are valid and the time index visual elements are relevant. Otherwise there is no time in the simulation and the elements are unnecessary. + * note it is always (as far as I can tell) the case that either all of the time indices are -1 (not being set by any object) or none are -1 (all are set by some object), hence + * we need only check the first time index */ + ifstream data(filename_tracks.toUtf8().constData(), ios::binary); + //define the number of bytes needed to skip to reach the first time index + int skipsize=sizeof(int)+sizeof(int)+sizeof(EGS_ParticleTrack::ParticleInfo); + //go to the position of the first time index and read it in + data.seekg(skipsize, ios::beg); + EGS_Float time; + data.read((char *)&time,sizeof(EGS_Float)); + data.close(); + //check that the time index is between 0 and 1. if so return true (thus setting hasdynamic true). Otherwise return false (similarly setting hasdynamic false) + if (time>=0 && time<1) { + return true; + } + else { + return false; + } +} + +void GeometryViewControl::timeObjectVisibility() { + //this function is used to make the time index objects visible or hidden as necessary. It will use the hasdynamic boolean, which indicates if they have any relevance in the uploaded files + if (!hasDynamic) { //if false all objects are hidden and the particle index spin boxes are allowed to be manually edited + button_timereset->hide(); + button_timeplay->hide(); + slider_timeindex->hide(); + spin_timewindow->hide(); + label_timewindow->hide(); + spin_timeindex->hide(); + label_timeindex->hide(); + groupBox_time->hide(); + spin_numTimeSteps->hide(); + label_numTimeSteps->hide(); + spin_tmaxe->setReadOnly(false); + spin_tmine->setReadOnly(false); + spin_tmaxpo->setReadOnly(false); + spin_tminpo->setReadOnly(false); + spin_tmaxp->setReadOnly(false); + spin_tminp->setReadOnly(false); + } + else { //if trueall objects are made visible and the particle index spin boxes cannot be manually edited + /* disable manual spinbox editing as the particle min and max indices will be determined by the time window and current time index (either from slider or spinbox). It is simpler in this case to simply + * have them modified only by interactions with time index objects to avoid any bugs or inconsistensies between the visual display and the dashboard settings*/ + spin_tmaxe->setReadOnly(true); + spin_tmine->setReadOnly(true); + spin_tmaxpo->setReadOnly(true); + spin_tminpo->setReadOnly(true); + spin_tmaxp->setReadOnly(true); + spin_tminp->setReadOnly(true); + button_timereset->show(); + button_timeplay->show(); + slider_timeindex->show(); + spin_timewindow->show(); + label_timewindow->show(); + spin_timeindex->show(); + label_timeindex->show(); + groupBox_time->show(); + spin_numTimeSteps->show(); + label_numTimeSteps->show(); + } + //has dynamic is true when a dynamic geometry is present, or when the tracks are being given some time index (exmaple due to a dynamic source or phasespace file) and include time = yes in input file (syncptracks file) + + +} + diff --git a/HEN_HOUSE/egs++/view/viewcontrol.h b/HEN_HOUSE/egs++/view/viewcontrol.h index 13c7d0e76..4959fd46d 100644 --- a/HEN_HOUSE/egs++/view/viewcontrol.h +++ b/HEN_HOUSE/egs++/view/viewcontrol.h @@ -26,6 +26,7 @@ # Contributors: Frederic Tessier # Ernesto Mainegra-Hing # Manuel Stoeckl +# Alexandre Demelo # ############################################################################### */ @@ -60,6 +61,7 @@ class GeometryViewControl : public QMainWindow, public Ui::GeometryViewControl { virtual void setFilename(QString str); virtual void setTracksFilename(QString str); + virtual void setTracksExtension(QString str); virtual void setCameraPosition(); virtual void setProjectionLineEdit(); virtual void setLightLineEdit(); @@ -76,6 +78,8 @@ class GeometryViewControl : public QMainWindow, public Ui::GeometryViewControl { virtual void updateRegionTable(int imed); virtual void updateAusgabObjects(bool loadUserDose=false); virtual void initColorSwatches(); + virtual bool hasValidTime(); + virtual void timeObjectVisibility(); public slots: @@ -142,7 +146,13 @@ public slots: virtual void changeTrackMaxP(int t); virtual void changeTrackMaxE(int t); virtual void changeTrackMaxPo(int t); - virtual void updateTracks(vector ntracks); + virtual void playTime(); + virtual void resetTime(); + virtual void slideTime(); + virtual void spinTime(); + virtual void updateNumTimeSteps(); + virtual void particleSlider(EGS_Float slidertime); + virtual void updateTracks(vector ntracks, vector timeindexlist_p, vector timeindexlist_e, vector timeindexlist_po); private: @@ -150,8 +160,13 @@ public slots: ImageWindow *gview; SaveImage *save_image; + vector timelist_p; //list of time indices of compressed particle tracks for photons + vector timelist_e; //list of time indices of compressed particle tracks for electrons + vector timelist_po; //list of time indices of compressed particle tracks for positrons + QString filename; QString filename_tracks; + QString tracks_extension; QString userDoseFile; int nmed; QRgb *m_colors; @@ -194,6 +209,8 @@ public slots: bool showPhotonTracks; bool showElectronTracks; bool showPositronTracks; + bool hasDynamic; //boolean to track whether or not to make time objects visible or hidden + bool isPlaying; vector show_regions; bool allowRegionSelection, energyScaling; diff --git a/HEN_HOUSE/egs++/view/viewcontrol.ui b/HEN_HOUSE/egs++/view/viewcontrol.ui index bec3933d2..7b9e93a4d 100644 --- a/HEN_HOUSE/egs++/view/viewcontrol.ui +++ b/HEN_HOUSE/egs++/view/viewcontrol.ui @@ -27,6 +27,7 @@ # # Contributors: Ernesto Mainegra-Hing # Frederic Tessier +# Alexandre Demelo # ############################################################################### @@ -36,8 +37,8 @@ 0 0 - 675 - 592 + 717 + 655 @@ -90,8 +91,8 @@ 0 0 - 671 - 563 + 713 + 626 @@ -130,7 +131,7 @@ - + QLayout::SetDefaultConstraint @@ -208,13 +209,13 @@ Particle tracks - - - - 1 + + + + Electrons - - 1 + + true @@ -228,6 +229,32 @@ + + + + 1 + + + 1 + + + + + + + + 0 + 0 + + + + 1 + + + 1 + + + @@ -238,18 +265,24 @@ - - - - Electrons + + + + 1 - - true + + 1 + + + 0 + 0 + + 1 @@ -258,8 +291,8 @@ - - + + 1 @@ -268,8 +301,14 @@ - - + + + + + 0 + 0 + + 1 @@ -278,29 +317,190 @@ - - - - 1 + + + + Time window + + + + + + + 3 + + + 2.000000000000000 + + + 0.001000000000000 - 1 + 0.010000000000000 - - - - 1 + + + + + 62 + 29 + - - 1 + + Reset + + + + + + + + 0 + 0 + + + + + 0 + 114 + + + + Time control + + + + + 10 + 30 + 231 + 21 + + + + + 0 + 0 + + + + 999 + + + 1 + + + Qt::Horizontal + + + + + + 180 + 50 + 62 + 26 + + + + + 62 + 26 + + + + + 62 + 26 + + + + Play + + + + + + 110 + 50 + 63 + 26 + + + + true + + + 3 + + + 0.999000000000000 + + + 0.001000000000000 + + + 0.000000000000000 + + + + + + 10 + 50 + 91 + 26 + + + + Time index + + + + + + 180 + 80 + 61 + 29 + + + + 1 + + + 1000000 + + + 1000 + + + 1000 + + + + + + 10 + 85 + 151 + 20 + + + + Number of time steps + + + + @@ -431,6 +631,14 @@ true + + + 0 + 0 + 232 + 68 + + @@ -1672,7 +1880,7 @@ p, li { white-space: pre-wrap; } 0 0 - 675 + 717 25 @@ -2654,5 +2862,101 @@ p, li { white-space: pre-wrap; } + + button_timeplay + pressed() + GeometryViewControl + playTime() + + + 488 + 321 + + + 358 + 327 + + + + + slider_timeindex + valueChanged(int) + GeometryViewControl + slideTime() + + + 608 + 321 + + + 358 + 327 + + + + + spin_timewindow + valueChanged(double) + GeometryViewControl + slideTime() + + + 590 + 241 + + + 358 + 327 + + + + + button_timereset + pressed() + GeometryViewControl + resetTime() + + + 658 + 241 + + + 358 + 327 + + + + + spin_timeindex + valueChanged(double) + GeometryViewControl + spinTime() + + + 594 + 334 + + + 358 + 327 + + + + + spin_numTimeSteps + valueChanged(int) + GeometryViewControl + updateNumTimeSteps() + + + 594 + 334 + + + 358 + 327 + + +