47 #include "../egs_gtransformed/egs_gtransformed.h"
63 const string EGS_AENVELOPE_LOCAL EGS_AEnvelope::allowed_base_geom_types[] = {
"EGS_cSpheres",
"EGS_cSphericalShell",
"EGS_XYZGeometry",
"EGS_RZ"};
65 static char EGS_AENVELOPE_LOCAL geom_class_msg[] =
"createGeometry(AEnvelope): %s\n";
66 static char EGS_AENVELOPE_LOCAL base_geom_keyword[] =
"base geometry";
67 static char EGS_AENVELOPE_LOCAL inscribed_geom_keyword[] =
"inscribed geometry";
68 static char EGS_AENVELOPE_LOCAL inscribed_geom_name_keyword[] =
"inscribed geometry name";
69 static char EGS_AENVELOPE_LOCAL transformations_keyword[] =
"transformations";
70 static char EGS_AENVELOPE_LOCAL type_keyword[] =
"type";
71 static char EGS_AENVELOPE_LOCAL transformation_keyword[] =
"transformation";
75 const vector<AEnvelopeAux> inscribed,
const string &Name,
bool debug,
string output_vc_file) :
76 EGS_BaseGeometry(Name), base_geom(base_geom), debug_info(debug), output_vc(output_vc_file) {
78 bool volcor_available = allowedBaseGeomType(base_geom->
getType());
79 bool volcor_requested = inscribed.size() > 0 && inscribed[0].vcopts->mode ==
CORRECT_VOLUME;
81 if (!volcor_available && volcor_requested) {
84 "EGS_AEnvelope:: Volume correction is not available for geometry type '%s (%s)'. "
85 "Geometry types must implement getVolume. Valid choices are:\n\t"
88 int end = (int)(
sizeof(allowed_base_geom_types)/
sizeof(string));
89 for (
int i=0; i < end; i++) {
90 msg += allowed_base_geom_types[i] +
" ";
92 msg +=
"\nPlease set `action = discover -or- correct and zero volume` or use a different base geometry type.\n";
98 nregbase = base_geom->
regions();
100 has_rho_scaling = getHasRhoScaling();
102 ninscribed = inscribed.size();
103 if (ninscribed == 0) {
104 egsFatal(
"EGS_AEnvelope: no inscribed geometries!\n");
112 for (
int geom_idx=0; geom_idx < ninscribed; geom_idx++) {
116 VCOptions *vcopt = inscribed[geom_idx].vcopts;
119 inscribed_geoms.push_back(transformed);
120 transforms.push_back(transform);
121 opts.push_back(vcopt);
123 int ninscribed_reg= transformed->
regions();
124 nreg += ninscribed_reg;
126 for (
int lreg = 0; lreg < ninscribed_reg; lreg++) {
128 local_to_global_reg[local] = global;
129 global_reg_to_local[global] = local;
134 geoms_in_region =
new vector<EGS_BaseGeometry *>[nregbase];
136 if (inscribed[0].vcopts->vc_file !=
"") {
137 vc_results =
loadFileResults(inscribed[0].vcopts, base_geom, inscribed_geoms, transforms);
138 egsInformation(
"loaded from %s\n", inscribed[0].vcopts->vc_file.c_str());
146 nreg_with_inscribed = 0;
147 for (volcor::RegionGeomSetT::iterator it=vc_results.regions_with_inscribed.begin();
148 it!=vc_results.regions_with_inscribed.end(); it++) {
149 nreg_with_inscribed++;
150 copy(it->second.begin(), it->second.end(), back_inserter(geoms_in_region[it->first]));
153 if (getNRegWithInscribed() == 0) {
154 egsFatal(
"EGS_AEnvelope:: Failed to find any regions with inscribed geometries\n");
157 if (output_vc==
"yes" || output_vc==
"text" || output_vc==
"gzip") {
158 writeVolumeCorrection();
168 EGS_AEnvelope::~EGS_AEnvelope() {
169 if (!base_geom->
deref()) {
177 int end = (int)(
sizeof(allowed_base_geom_types)/
sizeof(string));
179 for (
int i=0; i<end; i++) {
180 if (allowed_base_geom_types[i] == geom_type) {
188 int EGS_AEnvelope::getGlobalRegFromLocalReg(
EGS_BaseGeometry *g,
int local_reg) {
194 map<volcor::GeomRegPairT, int>::const_iterator glob_it = local_to_global_reg.find(local);
195 if (glob_it != local_to_global_reg.end()) {
196 return (*glob_it).second;
203 map<int, volcor::GeomRegPairT>::const_iterator glob_it = global_reg_to_local.find(ireg);
204 if (glob_it != global_reg_to_local.end()) {
205 return (*glob_it).second;
210 int EGS_AEnvelope::getNRegWithInscribed()
const {
211 return nreg_with_inscribed;
214 bool EGS_AEnvelope::isRealRegion(
int ireg)
const {
216 bool is_outside = ireg < 0 || ireg >=
nreg;
221 bool in_base_geom = ireg <
nregbase;
227 local = getLocalFromGlobalReg(ireg);
228 return local.first->isRealRegion(local.second);
233 bool EGS_AEnvelope::isInside(
const EGS_Vector &x) {
238 int EGS_AEnvelope::isWhere(
const EGS_Vector &x) {
240 int base_reg = base_geom->
isWhere(x);
243 vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(base_reg);
245 if (base_reg < 0 || geoms_in_reg.size()==0) {
251 for (vector<EGS_BaseGeometry *>::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
252 int inscribed_reg = (*geom)->isWhere(x);
253 if (inscribed_reg >= 0) {
254 return getGlobalRegFromLocalReg(*geom, inscribed_reg);
262 int EGS_AEnvelope::inside(
const EGS_Vector &x) {
266 int EGS_AEnvelope::medium(
int ireg)
const {
268 if (ireg < nregbase) {
269 return base_geom->
medium(ireg);
274 return local.first->medium(local.second);
277 vector<EGS_BaseGeometry *> EGS_AEnvelope::getGeomsInRegion(
int ireg) {
278 if (ireg < 0 || ireg >= nregbase) {
279 vector<EGS_BaseGeometry *> empty;
282 return geoms_in_region[ireg];
285 int EGS_AEnvelope::computeIntersections(
int ireg,
int n,
const EGS_Vector &X,
309 bool outside_envelope = ireg < 0;
310 bool in_base_geom = ireg >= 0 && ireg <
nregbase;
312 if (outside_envelope) {
316 ireg = howfar(ireg,x,u,t,&imed);
317 bool would_not_intersect = ireg < 0;
319 if (would_not_intersect) {
324 isections[0].
rhof = 1;
325 isections[0].
ireg = -1;
326 isections[0].
imed = -1;
351 isections[isec_idx].
imed = imed;
352 isections[isec_idx].
ireg = ireg;
353 isections[isec_idx].
rhof = getRelativeRho(ireg);
360 ireg = base_geom->
howfar(ireg,x,u,t,&imed);
363 vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(ireg);
364 for (vector<EGS_BaseGeometry *>::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
366 int inscribed_reg = (*geom)->howfar(-1, x, u, t, &imed);
368 bool hits_inscribed = inscribed_reg >= 0;
370 if (hits_inscribed) {
371 ireg = getGlobalRegFromLocalReg(*geom, inscribed_reg);
376 isections[isec_idx++].
t = ttot;
398 int max_isec = n - isec_idx;
400 int ninter_sec = local.first->computeIntersections(local.second, max_isec,
401 x, u, &isections[isec_idx]);
404 int nmax = ninter_sec >= 0 ? ninter_sec + isec_idx : n;
405 for (
int i=isec_idx; i < nmax; i++) {
406 isections[i].
ireg = getGlobalRegFromLocalReg(local.first, isections[i].
ireg);
407 isections[i].
t += ttot;
411 if (ninter_sec < 0) {
415 isec_idx += ninter_sec;
421 t = isections[isec_idx-1].
t - ttot;
423 ttot = isections[isec_idx-1].
t;
431 imed = base_geom->
medium(ireg);
445 bool in_base_geom = ireg <
nregbase;
449 else if (base_geom->
regions() == 1) {
453 int ir = base_geom->
isWhere(x);
461 EGS_Float &t,
int *newmed,
EGS_Vector *normal) {
463 bool inside_geom = ireg >= 0;
466 bool inside_base_geom = ireg <
nregbase;
468 if (inside_base_geom) {
470 int base_reg = base_geom->
howfar(ireg,x,u,t,newmed,normal);
471 vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(ireg);
473 bool no_inscribed_in_reg = geoms_in_reg.size() == 0;
474 if (no_inscribed_in_reg) {
478 int inscribed_global = -1;
479 for (vector<EGS_BaseGeometry *>::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
480 int local_reg = (*geom)->howfar(-1, x, u, t, newmed, normal);
481 bool hits_inscribed = local_reg >= 0;
482 if (hits_inscribed) {
483 inscribed_global = getGlobalRegFromLocalReg(*geom, local_reg);
487 bool hit_inscribed_first = inscribed_global >= 0;
489 if (hit_inscribed_first) {
490 return inscribed_global;
501 int new_local = local.first->howfar(local.second, x, u, t, newmed, normal);
502 bool hit_boundary_in_inscribed = new_local >=0;
503 if (hit_boundary_in_inscribed) {
504 return getGlobalRegFromLocalReg(local.first, new_local);
509 int inew = base_geom->
isWhere(x + u*t);
510 if (inew >= 0 && newmed) {
511 *newmed = base_geom->
medium(inew);
519 int new_base_reg = base_geom->
howfar(ireg,x,u,t,newmed,normal);
520 vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(new_base_reg);
521 bool has_geoms_in_reg = geoms_in_reg.size() > 0;
524 if (new_base_reg >= 0 && has_geoms_in_reg) {
527 for (vector<EGS_BaseGeometry *>::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
528 int new_local_reg = (*geom)->isWhere(x+u*t);
529 bool landed_in_inscribed = new_local_reg >= 0;
530 if (landed_in_inscribed) {
532 *newmed = (*geom)->medium(new_local_reg);
534 return getGlobalRegFromLocalReg(*geom, new_local_reg);
544 EGS_Float EGS_AEnvelope::hownear(
int ireg,
const EGS_Vector &x) {
545 bool inside_envelope = ireg >= 0;
547 if (!inside_envelope) {
548 return base_geom->
hownear(ireg, x);
551 bool in_base_geom = ireg <
nregbase ;
554 return local.first->hownear(local.second, x);
557 EGS_Float tmin = base_geom->
hownear(ireg, x);
561 vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(ireg);
562 for (vector<EGS_BaseGeometry *>::iterator geom=geoms_in_reg.begin(); geom!= geoms_in_reg.end(); ++geom) {
563 EGS_Float local_t = (*geom)->hownear(-1, x);
564 if (local_t < tmin) {
574 bool EGS_AEnvelope::hasBooleanProperty(
int ireg, EGS_BPType prop)
const {
575 if (ireg >= 0 && ireg <
nreg) {
576 if (ireg < nregbase) {
580 return local.first->hasBooleanProperty(local.second, prop);
585 void EGS_AEnvelope::setBooleanProperty(EGS_BPType prop) {
586 setPropertyError(
"setBooleanProperty()");
589 void EGS_AEnvelope::addBooleanProperty(
int bit) {
590 setPropertyError(
"addBooleanProperty()");
593 void EGS_AEnvelope::setBooleanProperty(EGS_BPType prop,
int start,
int end,
int step) {
594 setPropertyError(
"setBooleanProperty()");
597 void EGS_AEnvelope::addBooleanProperty(
int bit,
int start,
int end,
int step) {
598 setPropertyError(
"addBooleanProperty()");
601 int EGS_AEnvelope::getMaxStep()
const {
609 EGS_Float EGS_AEnvelope::getVolume(
int ireg) {
614 else if (ireg < nregbase) {
620 return local.first->getVolume(local.second);
623 EGS_Float EGS_AEnvelope::getCorrectionRatio(
int ireg) {
625 if (0 <= ireg && ireg < nregbase) {
637 void EGS_AEnvelope::printInfo()
const {
643 for (
int ir=0; ir <
nregbase; ir++) {
645 egsInformation(
" volume of region %d was corrected from %.5E g to %.5E g\n",
650 for (
int ir=0; ir <
nregbase; ir++) {
651 if (geoms_in_region[ir].size() > 0) {
652 egsInformation(
" region %d has %d inscribed geometries\n", ir, geoms_in_region[ir].size());
664 egsInformation(
"=======================================================\n");
667 void EGS_AEnvelope::writeVCToFile(ostream &out) {
669 vector<int> to_write;
671 for (
int i=0; i < base_geom->
regions(); i++) {
674 bool has_correction = fabs(cor-uncor) > 1E-8;
675 if (has_correction) {
676 to_write.push_back(i);
680 size_t nrecords = to_write.size();
681 out << nrecords << endl;
683 for (
size_t i=0; i<to_write.size(); i++) {
684 int ir = to_write[i];
688 vector<EGS_BaseGeometry *> geoms_in_reg = getGeomsInRegion(ir);
689 vector<int> geom_idxs_in_reg;
690 for (
size_t i=0; i<geoms_in_reg.size(); i++) {
693 geom_idxs_in_reg.push_back(pos);
697 int ninscribed = (int)geom_idxs_in_reg.size();
699 out << ir <<
" " << cor;
702 for (
size_t gidx = 0; gidx < geom_idxs_in_reg.size(); gidx++) {
703 out <<
" " <<geom_idxs_in_reg[gidx];
710 void EGS_AEnvelope::writeVolumeCorrection() {
712 egsInformation(
"\nVolume Correction Output File \n%s\n",
string(80,
'=').c_str());
714 string gname = base_geom->
getName();
715 string fname = gname+
".autoenv.volcor";
716 fname += (output_vc ==
"gzip" ?
".gz" :
"");
717 string mode = (output_vc ==
"gzip" ?
"gzip" :
"text");
720 "Writing volume correction file for %s to %s file %s\n",
721 gname.c_str(), mode.c_str(), fname.c_str());
723 if (output_vc ==
"gzip") {
725 ogzstream outgz(fname.c_str());
726 writeVCToFile(outgz);
731 ofstream out(fname.c_str());
740 bool EGS_AEnvelope::getHasRhoScaling() {
747 for (
size_t geom_idx=0; geom_idx <
inscribed_geoms.size(); geom_idx++) {
758 egsWarning(
"EGS_AEnvelope::setMedia: don't use this method. Use the\n"
759 " setMedia() methods of the geometry objects that make up this geometry\n");
763 void EGS_AEnvelope::setRelativeRho(
int start,
int end, EGS_Float rho) {
768 void EGS_AEnvelope::setRelativeRho(
EGS_Input *) {
769 egsWarning(
"EGS_AEnvelope::setRelativeRho(): don't use this method."
770 " Use the\n setRelativeRho methods of the geometry objects that make up"
775 EGS_Float EGS_AEnvelope::getRelativeRho(
int ireg)
const {
776 if (ireg < 0 || ireg >=
nreg) {
779 if (ireg < nregbase) {
785 return local.first->getRelativeRho(local.second);
796 const vector<AEnvelopeAux> inscribed,
const string &Name,
bool debug,
string output_vc_file):
797 EGS_AEnvelope(base_geom, inscribed, Name, debug, output_vc_file) {
809 vector<EGS_BaseGeometry *> EGS_ASwitchedEnvelope::getGeomsInRegion(
int ireg) {
811 vector<EGS_BaseGeometry *> active_in_reg;
812 if ((0 <= ireg) && (ireg < nregbase)) {
814 geoms_in_region[ireg].begin(), geoms_in_region[ireg].end(),
815 active_inscribed.begin(), active_inscribed.end(),
816 back_inserter(active_in_reg)
820 return active_in_reg;
826 active_inscribed.clear();
827 active_inscribed = geoms;
829 if (geoms.size() <= 0) {
844 active_inscribed.clear();
846 for (
size_t i = 0; i < geom_indexes.size(); i++) {
847 int idx = geom_indexes[i];
848 if (idx < 0 || idx >= ninscribed) {
849 egsFatal(
"EGS_ASwitchedEnvelope:: %d is not a valid geometry index\n", idx);
854 cur_ptr = geom_indexes[0];
861 if ((0 <= ireg) && (ireg < nregbase)) {
862 size_t nactive = getGeomsInRegion(ireg).size();
872 if ((0 <= ireg) && (ireg < nregbase)) {
873 size_t nactive = getGeomsInRegion(ireg).size();
874 size_t ntotal = EGS_AEnvelope::getGeomsInRegion(ireg).size();
875 return nactive < ntotal;
883 vector<EGS_BaseGeometry *> to_activate;
886 cur_ptr = inscribed_index;
894 vector<EGS_BaseGeometry *>::iterator loc = find(
895 active_inscribed.begin(), active_inscribed.end(),
898 if (loc != active_inscribed.end()) {
899 active_inscribed.erase(loc);
904 cur_ptr = cur_ptr == ninscribed -1 ? 0 : cur_ptr + 1;
905 vector<EGS_BaseGeometry *> to_activate;
920 egsWarning(
"Invalid transform input given\n");
923 transforms.push_back(transform);
944 vector<string> debug_choices;
945 debug_choices.push_back(
"no");
946 debug_choices.push_back(
"yes");
947 debug = input->
getInput(
"print debug info", debug_choices, 0);
949 int output_vc_file_choice;
950 vector<string> vc_file_choices;
951 vc_file_choices.push_back(
"no");
952 vc_file_choices.push_back(
"yes");
953 vc_file_choices.push_back(
"text");
954 vc_file_choices.push_back(
"gzip");
955 output_vc_file_choice = input->
getInput(
"output volume correction file", vc_file_choices, 0);
956 string output_vc_file = vc_file_choices[output_vc_file_choice];
959 if (output_vc_file ==
"gzip") {
961 "GZip file output requested but not compiled with gzstream.\n"
962 "Please recompile with gzstream support.\n"
970 string base_geom_name;
971 int err = input->
getInput(base_geom_keyword, base_geom_name);
973 egsWarning(geom_class_msg, (
"'"+
string(base_geom_keyword)+
"' input not found").c_str());
978 err = input->
getInput(type_keyword, type);
985 egsWarning(geom_class_msg, (
"Unable to find geometry '"+base_geom_name+
"'").c_str());
990 if (!inscribed_input) {
991 egsWarning(geom_class_msg, (
"Missing '"+
string(inscribed_geom_keyword)+
"' input item").c_str());
995 string inscribed_geom_name;
996 err = inscribed_input->
getInput(inscribed_geom_name_keyword, inscribed_geom_name);
998 egsWarning(geom_class_msg, (
"'"+
string(inscribed_geom_name_keyword)+
"' input not found").c_str());
1003 if (!inscribed_geom) {
1004 egsWarning(geom_class_msg, (
"Unable to find geometry '"+inscribed_geom_name+
"'").c_str());
1008 vector<EGS_AffineTransform *> transforms;
1018 if (!vcopts->
valid) {
1019 egsWarning(geom_class_msg,
"Missing or invalid 'region discovery' input item");
1023 delete inscribed_input;
1026 vector<AEnvelopeAux> inscribed;
1027 if (transforms.size()>0) {
1028 for (
size_t i=0; i < transforms.size(); i++) {
1029 inscribed.push_back(
AEnvelopeAux(inscribed_geom, transforms[i], vcopts));
1034 inscribed.push_back(
AEnvelopeAux(inscribed_geom, unityt, vcopts));
1038 if (type ==
"EGS_ASwitchedEnvelope") {
1042 result =
new EGS_AEnvelope(base_geom, inscribed,
"", (
bool)debug, output_vc_file);
virtual const string & getType() const =0
Get the geometry type.
virtual bool isInside(const EGS_Vector &x)=0
Is the position x inside the geometry?
virtual bool hasBooleanProperty(int ireg, EGS_BPType prop) const
Is the boolean property prop set for region ireg ?
A fast envelope geometry with automatic region detection.
int regions() const
Returns the number of local regions in this geometry.
int deref()
Decrease the reference count to this geometry.
An envelope geometry with automatic inscribed region detection (inspired by EGS_FastEnvelope) ...
A helper class for initializing auto envelopes.
virtual int medium(int ireg) const
Returns the medium index in region ireg.
virtual void printInfo() const
Print information about this geometry.
EGS_Float t
distance to next region boundary
This geometry type allows you to activate and deactivate inscribed geometries in custom egspp user co...
void outputResults() const
virtual EGS_Float howfarToOutside(int ireg, const EGS_Vector &x, const EGS_Vector &u)
EGS_Float rhof
relative mass density in that region
virtual int getMaxStep() const
Returns the maximum number of steps through the geometry.
A class representing 3D vectors.
pair< EGS_BaseGeometry *, int > GeomRegPairT
bool hasRhoScaling() const
Does this geometry object have a mass density scaling feature?
int nregbase
Number of regions in the base geometry.
virtual int isWhere(const EGS_Vector &x)=0
In which region is poisition x?
vector< EGS_BaseGeometry * > inscribed_geoms
The inscribed geometries.
void activateByIndex(int inscribed_index)
virtual int howfar(int ireg, const EGS_Vector &x, const EGS_Vector &u, EGS_Float &t, int *newmed=0, EGS_Vector *normal=0)=0
Calculate the distance to a boundary from x along the direction u.
EGS_GLIB_EXPORT EGS_BaseGeometry * createGeometry(EGS_Input *input)
Base geometry class. Every geometry class must be derived from EGS_BaseGeometry.
vector< EGS_AffineTransform * > transforms
The inscribed geometries.
bool hasActiveGeom(int ireg)
VCResults findRegionsWithInscribed(VCOptions *opts, EGS_BaseGeometry *base, vector< EGS_BaseGeometry * > inscribed, vector< EGS_AffineTransform * > transforms)
Run the MC simulation.
bool isConvex() const
Is the geometry convex?
void setActiveGeometries(vector< EGS_BaseGeometry * > geoms)
vector< EGS_Float > uncorrected_volumes
static EGS_BaseGeometry * getGeometry(const string &Name)
Get a pointer to the geometry named Name.
VCResults loadFileResults(VCOptions *opts, EGS_BaseGeometry *base, vector< EGS_BaseGeometry * > inscribed, vector< EGS_AffineTransform * > transforms)
Load volume corrections from external file.
EGS_InfoFunction EGS_EXPORT egsFatal
Always use this function for reporting fatal errors.
static bool allowedBaseGeomType(const string &geom_type)
function for checking whether a given geometry type is allowed to be used as a base geometry ...
EGS_InfoFunction EGS_EXPORT egsInformation
Always use this function for reporting the progress of a simulation and any other type of information...
EGS_BaseGeometry * base_geom
The envelope geometry.
void setName(EGS_Input *inp)
Set the name of the geometry from the input inp.
vector< EGS_Float > corrected_volumes
bool hasInactiveGeom(int ireg)
virtual EGS_Float hownear(int ireg, const EGS_Vector &x)=0
Calculate the distance to a boundary for position x in any direction.
virtual EGS_Float getRelativeRho(int ireg) const
Get the relative mass density in region ireg.
void deactivateByIndex(int inscribed_index)
bool has_rho_scaling
Does this geometry have relative mass density scvaling?
int ref()
Increase the reference count to this geometry.
virtual bool isRealRegion(int ireg) const
Returnes true if ireg is a real region, false otherwise.
static string type
Geometry type.
Volume correction initialization helper class.
const string & getName() const
Get the name of this geometry.
int nreg
Number of local regions in this geometry.
void setMedia(EGS_Input *, int, const int *)
Don't set media for an envelope geometry.
static vector< EGS_AffineTransform * > createTransforms(EGS_Input *inpt)
Take a block of transformations and return vector of EGS_AffineTransforms.
int ninscribed
Number of regions in the base geometry.
void setActiveByIndex(int inscribed_index)
static string type
Geometry type.
EGS_InfoFunction EGS_EXPORT egsWarning
Always use this function for reporting warnings.