libraries/HDGEOMETRY/DGeometry.cc

Bug Summary

File:	/home/sdobbs/work/clang/halld_recon/src/libraries/HDGEOMETRY/DGeometry.cc
Warning:	line 2093, column 24 The left operand of '+' is a garbage value

Annotated Source Code

Press '?' to see keyboard shortcuts

Show analyzer invocation

clang -cc1 -cc1 -triple x86_64-unknown-linux-gnu -analyze -disable-free -main-file-name DGeometry.cc -analyzer-store=region -analyzer-opt-analyze-nested-blocks -analyzer-checker=core -analyzer-checker=apiModeling -analyzer-checker=unix -analyzer-checker=deadcode -analyzer-checker=cplusplus -analyzer-checker=security.insecureAPI.UncheckedReturn -analyzer-checker=security.insecureAPI.getpw -analyzer-checker=security.insecureAPI.gets -analyzer-checker=security.insecureAPI.mktemp -analyzer-checker=security.insecureAPI.mkstemp -analyzer-checker=security.insecureAPI.vfork -analyzer-checker=nullability.NullPassedToNonnull -analyzer-checker=nullability.NullReturnedFromNonnull -analyzer-output plist -w -setup-static-analyzer -mrelocation-model pic -pic-level 2 -fhalf-no-semantic-interposition -mframe-pointer=none -fmath-errno -fno-rounding-math -mconstructor-aliases -munwind-tables -target-cpu x86-64 -tune-cpu generic -fno-split-dwarf-inlining -debugger-tuning=gdb -resource-dir /w/halld-scifs17exp/home/sdobbs/clang/llvm-project/install/lib/clang/12.0.0 -D HAVE_CCDB -D HAVE_RCDB -D HAVE_EVIO -D HAVE_TMVA=1 -D RCDB_MYSQL=1 -D RCDB_SQLITE=1 -D SQLITE_USE_LEGACY_STRUCT=ON -I .Linux_CentOS7.7-x86_64-gcc4.8.5/libraries/HDGEOMETRY -I libraries/HDGEOMETRY -I . -I libraries -I libraries/include -I /w/halld-scifs17exp/home/sdobbs/clang/halld_recon/Linux_CentOS7.7-x86_64-gcc4.8.5/include -I external/xstream/include -I /usr/include/tirpc -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/root/root-6.08.06/include -I /w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/ccdb/ccdb_1.06.06/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/rcdb/rcdb_0.06.00/cpp/include -I /usr/include/mysql -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/sqlitecpp/SQLiteCpp-2.2.0^bs130/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/sqlite/sqlite-3.13.0^bs130/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/hdds/hdds-4.9.0/Linux_CentOS7.7-x86_64-gcc4.8.5/src -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/xerces-c/xerces-c-3.1.4/include -I /group/halld/Software/builds/Linux_CentOS7.7-x86_64-gcc4.8.5/evio/evio-4.4.6/Linux-x86_64/include -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5 -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/x86_64-redhat-linux -internal-isystem /usr/lib/gcc/x86_64-redhat-linux/4.8.5/../../../../include/c++/4.8.5/backward -internal-isystem /usr/local/include -internal-isystem /w/halld-scifs17exp/home/sdobbs/clang/llvm-project/install/lib/clang/12.0.0/include -internal-externc-isystem /include -internal-externc-isystem /usr/include -O2 -std=c++11 -fdeprecated-macro -fdebug-compilation-dir /home/sdobbs/work/clang/halld_recon/src -ferror-limit 19 -fgnuc-version=4.2.1 -fcxx-exceptions -fexceptions -vectorize-loops -vectorize-slp -analyzer-output=html -faddrsig -o /tmp/scan-build-2021-01-21-110224-160369-1 -x c++ libraries/HDGEOMETRY/DGeometry.cc

libraries/HDGEOMETRY/DGeometry.cc

→

1// $Id$
2//
3//    File: DGeometry.cc
4// Created: Thu Apr  3 08:43:06 EDT 2008
5// Creator: davidl (on Darwin swire-d95.jlab.org 8.11.1 i386)
6//

8#include <algorithm>
9using namespace std;

11#include "DGeometry.h"
12#include <JANA/JGeometryXML.h>
13#include "FDC/DFDCWire.h"
14#include "FDC/DFDCGeometry.h"
15#include <ansi_escape((char)0x1b).h>

17using namespace std;

19#ifndef M_TWO_PI6.28318530717958647692
20#define M_TWO_PI6.28318530717958647692 6.28318530717958647692
21#endif

23//---------------------------------
24// DGeometry    (Constructor)
25//---------------------------------
26DGeometry::DGeometry(JGeometry *jgeom, DApplication *dapp, int32_t runnumber)
27{
this->jgeom = jgeom;
this->dapp = dapp;
this->bfield = NULL__null; // don't ask for B-field object before we're asked for it
this->runnumber = runnumber;
this->materialmaps_read = false;
this->materials_read = false;

pthread_mutex_init(&bfield_mutex, NULL__null);
pthread_mutex_init(&materialmap_mutex, NULL__null);
pthread_mutex_init(&materials_mutex, NULL__null);

ReadMaterialMaps();
40}

42//---------------------------------
43// ~DGeometry    (Destructor)
44//---------------------------------
45DGeometry::~DGeometry()
46{
pthread_mutex_lock(&materials_mutex);
for(unsigned int i=0; i<materials.size(); i++)delete materials[i];
materials.clear();
pthread_mutex_unlock(&materials_mutex);

pthread_mutex_lock(&materialmap_mutex);
for(unsigned int i=0; i<materialmaps.size(); i++)delete materialmaps[i];
materialmaps.clear();
pthread_mutex_unlock(&materialmap_mutex);
56}

58//---------------------------------
59// GetBfield
60//---------------------------------
61DMagneticFieldMap* DGeometry::GetBfield(void) const
62{
pthread_mutex_lock(&bfield_mutex);
if(bfield == NULL__null) bfield = dapp->GetBfield(runnumber);
pthread_mutex_unlock(&bfield_mutex);

return bfield;
68}

70//---------------------------------
71// GetLorentzDeflections
72//---------------------------------
73DLorentzDeflections* DGeometry::GetLorentzDeflections(void)
74{
return dapp->GetLorentzDeflections(runnumber);
76}

78//---------------------------------
79// GetMaterialMapVector
80//---------------------------------
81vector<DMaterialMap*> DGeometry::GetMaterialMapVector(void) const
82{
83//	ReadMaterialMaps();

return materialmaps;
86}

88//---------------------------------
89// ReadMaterialMaps
90//---------------------------------
91void DGeometry::ReadMaterialMaps(void) const
92{
/// This gets called by several "FindMat" methods below. It
/// will return immediately if the material maps have already
/// been read in. If they have not been read in, then it reads
/// them and sets a flag so that they are only read in once.
///
/// Orginally, this code resided in the constructor, but was
/// moved here so that programs not requiring the material
/// maps could start up quicker and skip reading them in altogether.

// Lock mutex so we can check flag to see if maps have already
// been read in
pthread_mutex_lock(&materialmap_mutex);
if(materialmaps_read){
// Maps are already read. Unlock mutex and return now
pthread_mutex_unlock(&materialmap_mutex);
return;
}

JCalibration * jcalib = dapp->GetJCalibration(runnumber);
if(!jcalib){
_DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<113<<" "<<"ERROR:  Unable to get JCalibration object!"<<endl;
pthread_mutex_unlock(&materialmap_mutex);
return;
}

// Get a list of all namepaths so we can find all of the material maps.
// We want to read in all material maps with a standard naming scheme
// so the number and coverage of the piece-wise maps can be changed
// without requiring recompilation.
//vector<DMaterialMap*> material_maps;
vector<string> namepaths;
jcalib->GetListOfNamepaths(namepaths);
vector<string> material_namepaths;
for(unsigned int i=0; i<namepaths.size(); i++){
if(namepaths[i].find("Material/material_map")==0)material_namepaths.push_back(namepaths[i]);
}

// Sort alphabetically so user controls search sequence by map name
sort(material_namepaths.begin(), material_namepaths.end());

// Inform user what's happening
if(material_namepaths.size()==0){
jerr<<"No material maps found in calibration DB!!"<<endl;
pthread_mutex_unlock(&materialmap_mutex);
return;
}
jout<<"Found "<<material_namepaths.size()<<" material maps in calib. DB"<<endl;

if(false){ // save this to work off configuration parameter
jout<<"Will read in the following:"<<endl;
for(unsigned int i=0; i<material_namepaths.size(); i++){
	jout<<"  "<<material_namepaths[i]<<endl;
}
}

// Actually read in the maps
uint32_t Npoints_total=0;
//cout<<endl; // make empty line so material map can overwrite it below
for(unsigned int i=0; i<material_namepaths.size(); i++){
// DMaterialMap constructor prints line so we conserve real
// estate by having each recycle the line
//cout<<ansi_up(1)<<string(85, ' ')<<"\r";
DMaterialMap *mat = new DMaterialMap(material_namepaths[i], jcalib);
if( ! mat->IS_VALID ) {
	// This particular map may not exist for this run/variation
	// (e.g. CPP uses maps downstream of TOF)
	delete mat;
	continue;
}
materialmaps.push_back(mat);
Npoints_total += (unsigned int)(mat->GetNr()*mat->GetNz());
}
//cout<<ansi_up(1)<<string(85, ' ')<<"\r";
jout<<"Read in "<<materialmaps.size()<<" material maps for run "<<runnumber<<" containing "<<Npoints_total<<" grid points total"<<endl;

// Set flag that maps have been read and unlock mutex
materialmaps_read = true;
pthread_mutex_unlock(&materialmap_mutex);
171}

173//---------------------------------
174// FindNodes
175//---------------------------------
176void DGeometry::FindNodes(string xpath, vector<xpathparsed_t> &matched_xpaths) const
177{
/// Find all nodes that match the specified xpath and return them as
/// fully parsed lists of the nodes and attributes.
///
/// The matched_xpaths variable has 4 levels of STL containers nested
/// together! The node_t data type contains 2 of them and represents
/// a single tag with the "first" member being the tag name
/// and the "second" member being a map of all of the attributes
/// of that tag along with their values.
///
/// The xpathparsed_t data type is a STL vector of node_t objects
/// that comprises a complete xpath. The data type of matched_xpaths
/// is a STL vector if xpathparsed_t objects and so comprises a
/// list of complete xpaths.

/// We do this by calling the GetXPaths() method of JGeometry to get
/// a list of all xpaths. Then we pass all of those in to 
/// JGeometryXML::ParseXPath() to get a parsed list for each. This
/// is compared to the parsed values of the xpath passed to us
/// (also parsed by JGeometryXML::ParseXPath()) to find matches.

// Make sure matched_xpaths is empty
matched_xpaths.clear();

// Cast JGeometry into a JGeometryXML
JGeometryXML *jgeomxml = dynamic_cast<JGeometryXML*>(jgeom);

// Parse our target xpath
xpathparsed_t target;
string unused_string;
unsigned int unused_int;
jgeomxml->ParseXPath(xpath, target, unused_string, unused_int);

// Get all xpaths for current geometry source
vector<string> allxpaths;
jgeom->GetXPaths(allxpaths, JGeometry::attr_level_all);

// Loop over xpaths
for(unsigned int i=0; i<allxpaths.size(); i++);
216}

218//---------------------------------
219// FindMatALT1 - parameters for alt1 fitter.
220//---------------------------------
221jerror_t DGeometry::FindMatALT1(DVector3 &pos, DVector3 &mom,double &KrhoZ_overA, 
		double &rhoZ_overA,double &LnI,
		double &X0, double *s_to_boundary) const
224{
225//	ReadMaterialMaps();

for(unsigned int i=0; i<materialmaps.size(); i++){
jerror_t err = materialmaps[i]->FindMatALT1(pos,KrhoZ_overA, rhoZ_overA,LnI,X0);
if(err==NOERROR){
	// We found the material map containing this point. If a non-NULL 
	// pointer was passed in for s_to_boundary, then search through all
	// material maps above and including this one to find the estimated
	// distance to the nearest boundary the swimmer should step to. Maps
	// after this one would only be considered outside of this one so
	// there is no need to check them.
	if(s_to_boundary==NULL__null)return NOERROR;	// User doesn't want distance to boundary
	*s_to_boundary = 1.0E6;
	for(unsigned int j=0; j<=i; j++){
		double s = materialmaps[j]->EstimatedDistanceToBoundary(pos, mom);
		if(s<*s_to_boundary)*s_to_boundary = s;
	}
	return NOERROR;
}
}
return RESOURCE_UNAVAILABLE;
246}




251//---------------------------------
252// FindMatKalman - Kalman filter needs slightly different set of parms.
253//---------------------------------
254jerror_t DGeometry::FindMatKalman(const DVector3 &pos,const DVector3 &mom,
		  double &KrhoZ_overA, 
		  double &rhoZ_overA, 
		  double &LnI,double &Z,
		  double &chi2c_factor,double &chi2a_factor,
		  double &chi2a_corr,
		  unsigned int &last_index,
		  double *s_to_boundary) const
262{
263//	ReadMaterialMaps();

//last_index=0;
for(unsigned int i=last_index; i<materialmaps.size(); i++){
  jerror_t err = materialmaps[i]->FindMatKalman(pos,KrhoZ_overA,
				  rhoZ_overA,LnI,chi2c_factor,
				  chi2a_factor,chi2a_corr,Z);
  if(err==NOERROR){
    if(i==materialmaps.size()-1) last_index=0;
    else last_index=i;
    if(s_to_boundary==NULL__null)return NOERROR;	// User doesn't want distance to boundary

    *s_to_boundary = 1.0E6;
    // If we are in the main mother volume, search through all the maps for
    // the nearest boundary
    if(last_index==0){
for(unsigned int j=0; j<materialmaps.size();j++){
 double s = materialmaps[j]->EstimatedDistanceToBoundary(pos, mom);
 if(s<*s_to_boundary){
   *s_to_boundary = s;
 }
}
    }
    else{
// otherwise, we found the material map containing this point. 
double s = materialmaps[last_index]->EstimatedDistanceToBoundary(pos, mom);
if(s<*s_to_boundary)*s_to_boundary = s;
    }
    return NOERROR;
  }
}
     
295return RESOURCE_UNAVAILABLE;
296}

298jerror_t DGeometry::FindMatKalman(const DVector3 &pos,const DVector3 &mom,
		  double &KrhoZ_overA, 
		  double &rhoZ_overA, 
		  double &LnI,double &Z,
		  unsigned int &last_index,
		  double *s_to_boundary) const
304{
305//	ReadMaterialMaps();

//last_index=0;
for(unsigned int i=last_index; i<materialmaps.size(); i++){
  jerror_t err = materialmaps[i]->FindMatKalman(pos,KrhoZ_overA,
				  rhoZ_overA,LnI,Z);
  if(err==NOERROR){
    if(i==materialmaps.size()-1) last_index=0;
    else last_index=i;
    if(s_to_boundary==NULL__null)return NOERROR;	// User doesn't want distance to boundary

    *s_to_boundary = 1.0E6;
    // If we are in the main mother volume, search through all the maps for
    // the nearest boundary
    if(last_index==0){
for(unsigned int j=0; j<materialmaps.size();j++){
 double s = materialmaps[j]->EstimatedDistanceToBoundary(pos, mom);
 if(s<*s_to_boundary){
   *s_to_boundary = s;
 }
}
    }
    else{
// otherwise, we found the material map containing this point. 
double s = materialmaps[last_index]->EstimatedDistanceToBoundary(pos, mom);
if(s<*s_to_boundary)*s_to_boundary = s;
    }
    return NOERROR;
  }
}
     
336return RESOURCE_UNAVAILABLE;
337}

339//---------------------------------
340jerror_t DGeometry::FindMatKalman(const DVector3 &pos,
		  double &KrhoZ_overA, 
		  double &rhoZ_overA, 
		  double &LnI, double &Z,double &chi2c_factor, 
		  double &chi2a_factor, double &chi2a_corr,
		  unsigned int &last_index) const
346{
347//	ReadMaterialMaps();

//last_index=0;
for(unsigned int i=last_index; i<materialmaps.size(); i++){
  jerror_t err = materialmaps[i]->FindMatKalman(pos,KrhoZ_overA,
				  rhoZ_overA,LnI,
				  chi2c_factor,chi2a_factor,
				  chi2a_corr,Z);
  if(err==NOERROR){
    if(i==materialmaps.size()-1) last_index=0;
    else last_index=i;
    return err;
  }
}
     
return RESOURCE_UNAVAILABLE;
363}

365//---------------------------------
366// Get material properties needed for dEdx
367jerror_t DGeometry::FindMatKalman(const DVector3 &pos,
		  double &KrhoZ_overA, 
		  double &rhoZ_overA, 
		  double &LnI, double &Z,
		  unsigned int &last_index) const
372{
373//	ReadMaterialMaps();

//last_index=0;
for(unsigned int i=last_index; i<materialmaps.size(); i++){
  jerror_t err = materialmaps[i]->FindMatKalman(pos,KrhoZ_overA,
				  rhoZ_overA,LnI,Z);
  if(err==NOERROR){
    if(i==materialmaps.size()-1) last_index=0;
    else last_index=i;
    return err;
  }
}
     
return RESOURCE_UNAVAILABLE;
387}




392//---------------------------------
393// FindMat
394//---------------------------------
395jerror_t DGeometry::FindMat(DVector3 &pos, double &rhoZ_overA, double &rhoZ_overA_logI, double &RadLen) const
396{
397//	ReadMaterialMaps();

for(unsigned int i=0; i<materialmaps.size(); i++){
jerror_t err = materialmaps[i]->FindMat(pos, rhoZ_overA, rhoZ_overA_logI, RadLen);
if(err==NOERROR)return NOERROR;
}
return RESOURCE_UNAVAILABLE;
404}

406//---------------------------------
407// FindMat
408//---------------------------------
409jerror_t DGeometry::FindMat(DVector3 &pos, double &density, double &A, double &Z, double &RadLen) const
410{
411//	ReadMaterialMaps();

for(unsigned int i=0; i<materialmaps.size(); i++){
jerror_t err = materialmaps[i]->FindMat(pos, density, A, Z, RadLen);
if(err==NOERROR)return NOERROR;
}
return RESOURCE_UNAVAILABLE;
418}

420//---------------------------------
421// FindMatNode
422//---------------------------------
423//const DMaterialMap::MaterialNode* DGeometry::FindMatNode(DVector3 &pos) const
424//{
425//
426//}

428//---------------------------------
429// FindDMaterialMap
430//---------------------------------
431const DMaterialMap* DGeometry::FindDMaterialMap(DVector3 &pos) const
432{
433//	ReadMaterialMaps();

for(unsigned int i=0; i<materialmaps.size(); i++){
const DMaterialMap* map = materialmaps[i];
if(map->IsInMap(pos))return map;
}
return NULL__null;
440}

442//====================================================================
443// Convenience Methods
444//
445// Below are defined some methods to make it easy to extract certain
446// key values about the GlueX detector geometry from the XML source.
447// Note that one can still use the generic Get(string xpath, ...) 
448// methods. This just packages some of them up for convenience.
449//
450// The one real gotcha here is that these methods must be kept
451// in sync with the XML structure by hand. If volumes are renamed
452// or their location within the hierachy is modified, then these
453// routines will need to be modified as well. That, or course, is
454// also true if you are using the generic Get(...) methods.
455//
456// What these methods are useful for is when minor changes are made
457// to the XML (such as the locations of the FDC packages) they
458// are automatically reflected here.
459//====================================================================

461//---------------------------------
462// GetDMaterial
463//---------------------------------
464const DMaterial* DGeometry::GetDMaterial(string name) const
465{
/// Get a pointer to the DMaterial object with the specified name.
// Only fill materials member when one is actually requested
// and then, only fill it once.

// Lock mutex so we can check flag to see if maps have already
// been read in
pthread_mutex_lock(&materials_mutex);
if(!materials_read) GetMaterials();
pthread_mutex_unlock(&materials_mutex);

for(unsigned int i=0; i<materials.size(); i++){
if(materials[i]->GetName() == name)return materials[i];
}

//_DBG_<<"No material \""<<name<<"\" found ("<<materials.size()<<" materials defined)"<<endl;

return NULL__null;
483}

485//---------------------------------
486// GetMaterials
487//---------------------------------
488void DGeometry::GetMaterials(void) const
489{
/// Read in all of the materials from the geometry source and create
/// a DMaterial object for each one.

//=========== elements ===========
string filter = "//materials/element/real[@name=\"radlen\"]";

// Get list of all xpaths
vector<string> xpaths;
jgeom->GetXPaths(xpaths, JGeometry::attr_level_all, filter);

// Look for xpaths that have "/materials[" in them
for(unsigned int i=0; i<xpaths.size(); i++){
// Get start of "element" section
string::size_type pos = xpaths[i].find("/element[");
if(pos == string::npos)continue;

// Get name attribute
string::size_type start_pos = xpaths[i].find("@name=", pos);
start_pos = xpaths[i].find("'", start_pos);
string::size_type end_pos = xpaths[i].find("'", start_pos+1);
if(end_pos==string::npos)continue;
string name = xpaths[i].substr(start_pos+1, end_pos-(start_pos+1));

// Get values
char xpath[256];

double A;
sprintf(xpath,"//materials/element[@name='%s']/[@a]", name.c_str());
if(!Get(xpath, A))continue;

double Z;
sprintf(xpath,"//materials/element[@name='%s']/[@z]", name.c_str());
if(!Get(xpath, Z))continue;

double density;
sprintf(xpath,"//materials/element[@name='%s']/real[@name='density']/[@value]", name.c_str());
if(!Get(xpath, density))continue;

double radlen;
sprintf(xpath,"//materials/element[@name='%s']/real[@name='radlen']/[@value]", name.c_str());
if(!Get(xpath, radlen))continue;

DMaterial *mat = new DMaterial(name, A, Z, density, radlen);
materials.push_back(mat);

cout<<mat->toString();
}

//=========== composites ===========
filter = "//materials/composite[@name]";

// Get list of all xpaths
jgeom->GetXPaths(xpaths, JGeometry::attr_level_all, filter);

// Look for xpaths that have "/materials[" in them
for(unsigned int i=0; i<xpaths.size(); i++){
// Get start of "composite" section
string::size_type pos = xpaths[i].find("/composite[");
if(pos == string::npos)continue;

// Get name attribute
string::size_type start_pos = xpaths[i].find("@name=", pos);
start_pos = xpaths[i].find("'", start_pos);
string::size_type end_pos = xpaths[i].find("'", start_pos+1);
if(end_pos==string::npos)continue;
string name = xpaths[i].substr(start_pos+1, end_pos-(start_pos+1));

if(GetDMaterial(name))continue; // skip duplicates

// Get values
char xpath[256];

// We should calculate an effective A and Z .... but we don't
bool found_all=true;
double A=0;
double Z=0;

double density;
sprintf(xpath,"//materials/composite[@name='%s']/real[@name='density']/[@value]", name.c_str());
found_all &= Get(xpath, density);

double radlen;
sprintf(xpath,"//materials/composite[@name='%s']/real[@name='radlen']/[@value]", name.c_str());
found_all &= Get(xpath, radlen);

// If we didn't find the info we need (radlen and density) in the 
// composite tag itself. Try calculating them from the components
if(!found_all)found_all = GetCompositeMaterial(name, density, radlen);

// If we weren't able to get the values needed to make the DMaterial object
// then skip this one.
if(!found_all)continue;

DMaterial *mat = new DMaterial(name, A, Z, density, radlen);
materials.push_back(mat);

cout<<mat->toString();
}

materials_read = true;
590}

592//---------------------------------
593// GetCompositeMaterial
594//---------------------------------
595bool DGeometry::GetCompositeMaterial(const string &name, double &density, double &radlen) const
596{
// Get list of all xpaths with "addmaterial" and "fractionmass"
char filter[512];
sprintf(filter,"//materials/composite[@name='%s']/addmaterial/fractionmass[@fraction]", name.c_str());
vector<string> xpaths;
jgeom->GetXPaths(xpaths, JGeometry::attr_level_all, filter);

// Loop over components of this composite
604_DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<604<<" "<<"Components for compsite "<<name<<endl;
for(unsigned int i=0; i<xpaths.size(); i++){
// Get component material name
string::size_type start_pos = xpaths[i].find("@material=", 0);
start_pos = xpaths[i].find("'", start_pos);
string::size_type end_pos = xpaths[i].find("'", start_pos+1);
if(end_pos==string::npos)continue;
string mat_name = xpaths[i].substr(start_pos+1, end_pos-(start_pos+1));

// Get component mass fraction
start_pos = xpaths[i].find("fractionmass[", 0);
start_pos = xpaths[i].find("@fraction=", start_pos);
start_pos = xpaths[i].find("'", start_pos);
end_pos = xpaths[i].find("'", start_pos+1);
if(end_pos==string::npos)continue;
string mat_frac_str = xpaths[i].substr(start_pos+1, end_pos-(start_pos+1));
double fractionmass = atof(mat_frac_str.c_str());

_DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<622<<" "<<"   "<<xpaths[i]<<"  fractionmass="<<fractionmass<<endl;
}

return true;
626}

628//---------------------------------
629// GetTraversedMaterialsZ
630//---------------------------------
631//void DGeometry::GetTraversedMaterialsZ(double q, const DVector3 &pos, const DVector3 &mom, double z_end, vector<DMaterialStep> &materialsteps)
632//{
/// Find the materials traversed by a particle swimming from a specific
/// position with a specific momentum through the magnetic field until
/// it reaches a specific z-location. Energy loss is not included in
/// the swimming since this method itself is assumed to be one of the
/// primary means of determining energy loss. As such, one should not
/// pass in a value of z_end that is far from pos.
///
/// The vector materialsteps will be filled with DMaterialStep objects
/// corresponding to each of the materials traversed.
///
/// The real work here is done by the DMaterialStepper class

645//}

647//---------------------------------
648// GetCDCStereoWires
649//---------------------------------
650/// Extract the stereo wire data from the XML
651bool DGeometry::GetCDCStereoWires(unsigned int ring,unsigned int ncopy,
		  double zcenter,double dz,  
		  vector<vector<cdc_offset_t> >&cdc_offsets,
		  vector<DCDCWire*> &stereowires,
		  vector<double>&rot_angles,double dx,
		  double dy) const{
stringstream r_z_s,phi0_s,rot_s;

// Create search strings for the straw geometrical properties 
r_z_s << "//mposPhi[@volume='CDCstrawLong']/@R_Z/ring[@value='" << ring << "']";
phi0_s << "//mposPhi[@volume='CDCstrawLong']/@Phi0/ring[@value='" << ring << "']";
rot_s << "//mposPhi[@volume='CDCstrawLong']/@rot/ring[@value='" << ring << "']";

vector<double>r_z;
double phi0;
vector<double>rot;

// Extract data from the XML
if(!Get(r_z_s.str(), r_z)) return false; 
if(!Get(phi0_s.str(), phi0)) return false; 
if(!Get(rot_s.str(), rot)) return false; 

// Angular quantities
const double deg2rad=M_PI3.14159265358979323846/180.;
double dphi=2*M_PI3.14159265358979323846/double(ncopy);
phi0*=deg2rad;

// Extract data from close-packed straws
double stereo=0.,stereo_sign=1.;
stereo=deg2rad*rot[0];
if (stereo<0.) stereo_sign=-1.;

// Loop over the number of straws
for (unsigned int i=0;i<ncopy;i++){
   DCDCWire *w=new DCDCWire;
   double phi=phi0+double(i)*dphi;
   w->ring=ring;
   w->straw=i+1;

   // Find the nominal wire position and direction from the XML
   DVector3 origin,udir;
   origin.SetX(r_z[0]*cos(phi)+dx);
   origin.SetY(r_z[0]*sin(phi)+dy);    
   origin.SetZ(zcenter);

   // Here, we need to define a coordinate system for the wire
   // in which the wire runs along one axis. We call the directions
   // of the axes in this coordinate system s,t, and u with
   // the wire running in the "u" direction. The "s" direction
   // will be defined by the direction pointing from the beamline
   // to the midpoint of the wire.
   udir.SetXYZ(0.0, 0.0,1.0);	
   udir.RotateX(stereo);	
   udir.RotateZ(phi);     

   // Apply offsets in x and y
   double half_dz=0.5*dz;
   double x0=origin.x(),y0=origin.y();
   double ux=udir.x()/udir.z();
   double uy=udir.y()/udir.z();
   unsigned int ringid=ring-1;
   DVector3 downstream(x0+half_dz*ux+cdc_offsets[ringid][i].dx_d,
         y0+half_dz*uy+cdc_offsets[ringid][i].dy_d,
         zcenter+half_dz);
   DVector3 upstream(x0-half_dz*ux+cdc_offsets[ringid][i].dx_u,
         y0-half_dz*uy+cdc_offsets[ringid][i].dy_u,
         zcenter-half_dz);
   w->origin=0.5*(upstream+downstream);
   w->origin.RotateX(rot_angles[0]);
   w->origin.RotateY(rot_angles[1]);
   w->origin.RotateZ(rot_angles[2]);

   w->phi=w->origin.Phi();

   // Wire direction
   w->udir=downstream-upstream;
   w->udir.RotateX(rot_angles[0]);
   w->udir.RotateY(rot_angles[1]);
   w->udir.RotateZ(rot_angles[2]);

   // For derivatives
   w->udir_mag=w->udir.Mag();

   w->udir.SetMag(1.);
   w->stereo=stereo_sign*w->udir.Angle(DVector3(0.,0.,1.));
   // other directions for our wire coordinate system
   w->sdir=w->origin;
   w->sdir.SetMag(1.);
   w->tdir = w->udir.Cross(w->sdir);

   // Some values needed for alignment derivatives
   w->x0=dx; w->y0=dy; w->z0=zcenter;
   w->phiX=rot_angles[0]; w->phiY=rot_angles[1]; w->phiZ=rot_angles[2];
   w->r0=r_z[0]; w->phiStraw=phi; w->stereo_raw=stereo;

   stereowires.push_back(w);
}

return true;
750}

752//---------------------------------
753// GetCDCAxialWires
754//---------------------------------
755/// Extract the axial wire data from the XML
756bool DGeometry::GetCDCAxialWires(unsigned int ring,unsigned int ncopy,
    double zcenter,double dz,
    vector<vector<cdc_offset_t> >&cdc_offsets,
    vector<DCDCWire*> &axialwires,
    vector<double>&rot_angles,double dx,
    double dy) const{
 stringstream phi0_s,r_z_s;

 // Create search strings for the number of straws and the straw geometrical properties 
 phi0_s << "//mposPhi[@volume='CDCstrawShort']/@Phi0/ring[@value='" << ring << "']";
 r_z_s << "//mposPhi[@volume='CDCstrawShort']/@R_Z/ring[@value='" << ring << "']";

 double phi0;
 vector<double>r_z;

 // Extract the data from the XML
 if(!Get(phi0_s.str(), phi0)) return false; 
 if(!Get(r_z_s.str(), r_z)) return false;

 // Angular quantities
 double dphi=2*M_PI3.14159265358979323846/double(ncopy);
 phi0*=M_PI3.14159265358979323846/180.;

 // Loop over the number of straws
 for (unsigned int i=0;i<ncopy;i++){
    DCDCWire *w=new DCDCWire;
    double phi=phi0+double(i)*dphi;
    w->ring=ring;
    w->straw=i+1;

    // Find the nominal wire position from the XML
    double x0=r_z[0]*cos(phi)+dx;
    double y0=r_z[0]*sin(phi)+dy;

    // Apply offsets in x and y
    double half_dz=0.5*dz;
    unsigned int ringid=ring-1;
    DVector3 downstream(x0+cdc_offsets[ringid][i].dx_d,
          y0+cdc_offsets[ringid][i].dy_d,
          zcenter+half_dz);
    DVector3 upstream(x0+cdc_offsets[ringid][i].dx_u,
          y0+cdc_offsets[ringid][i].dy_u,
          zcenter-half_dz);
    w->origin=0.5*(upstream+downstream);
    w->origin.RotateX(rot_angles[0]);
    w->origin.RotateY(rot_angles[1]);
    w->origin.RotateZ(rot_angles[2]);
    w->phi=w->origin.Phi();

    // Here, we need to define a coordinate system for the wire
    // in which the wire runs along one axis. We call the directions
    // of the axes in this coordinate system s,t, and u with
    // the wire running in the "u" direction. The "s" direction
    // will be defined by the direction pointing from the beamline
    // to the midpoint of the wire.
    w->udir=downstream-upstream;
    w->udir.RotateX(rot_angles[0]);
    w->udir.RotateY(rot_angles[1]);
    w->udir.RotateZ(rot_angles[2]);

    // For derivatives
    w->udir_mag=w->udir.Mag();

    w->udir.SetMag(1.);

    w->stereo=w->udir.Angle(DVector3(0.,0.,1.));
    // other directions for our wire coordinate system
    w->sdir=w->origin;
    w->sdir.SetMag(1.);
    w->tdir = w->udir.Cross(w->sdir);

    // Some values needed for alignment derivatives
    w->x0=dx; w->y0=dy; w->z0=zcenter;
    w->phiX=rot_angles[0]; w->phiY=rot_angles[1]; w->phiZ=rot_angles[2];
    w->r0=r_z[0]; w->phiStraw=phi; w->stereo_raw=0.0;

    axialwires.push_back(w);
 }

 return true;
836}

838//---------------------------------
839// GetCDCWires
840//---------------------------------
841bool DGeometry::GetCDCWires(vector<vector<DCDCWire *> >&cdcwires) const{
 // Get nominal geometry from XML
 vector<double>cdc_origin;
 vector<double>cdc_length;
 Get("//posXYZ[@volume='CentralDC']/@X_Y_Z",cdc_origin);
 Get("//tubs[@name='STRW']/@Rio_Z",cdc_length);

 // Get CDC rotation angles
 vector<double>rot_angles;
 Get("//posXYZ[@volume='CentralDC']/@rot", rot_angles);
 rot_angles[0]*=M_PI3.14159265358979323846/180.;
 rot_angles[1]*=M_PI3.14159265358979323846/180.;
 rot_angles[2]*=M_PI3.14159265358979323846/180.; 

 double dX=0.0, dY=0.0, dZ=0.0;
 double dPhiX=0.0,dPhiY=0.0,dPhiZ=0.0;

 JCalibration * jcalib = dapp->GetJCalibration(runnumber);
 vector<map<string,double> >vals;
 if (jcalib->Get("CDC/global_alignment",vals)==false){
    map<string,double> &row = vals[0];
    dX=row["dX"];
    dY=row["dY"];
    dZ=row["dZ"];
    dPhiX=row["dPhiX"];
    dPhiY=row["dPhiY"];
    dPhiZ=row["dPhiZ"];
 }


 // Shift the CDC origin according to alignment
 cdc_origin[0]+=dX;
 cdc_origin[1]+=dY;
 cdc_origin[2]+=dZ;

 // Shift the CDC rotation according to the alignment
 rot_angles[0]+=dPhiX;
 rot_angles[1]+=dPhiY;
 rot_angles[2]+=dPhiZ;

 double zmin=cdc_origin[2];
 double zmax=zmin+cdc_length[2];
 double zcenter=0.5*(zmin+zmax);
 double L=zmax-zmin;

 // Get number of straws for each layer from the XML
 unsigned int numstraws[28];  
 stringstream ncopy_s;

 // Get number of straws for each ring
 for (unsigned int ring=1;ring<=28;ring++){
    // Create search string for the number of straws 
    ncopy_s << "//section[@name='CentralDC']/composition/mposPhi/@ncopy/ring[@value='" << ring << "']";
    Get(ncopy_s.str(),numstraws[ring-1]);
    ncopy_s.str("");
    ncopy_s.clear();
 }

 // Get straw-by-straw offsets from calibration database
 vector<cdc_offset_t>tempvec;
 vector<vector<cdc_offset_t> >cdc_offsets;

 if (jcalib->Get("CDC/wire_alignment",vals)==false){
    unsigned int straw_count=0,ring_count=0;
    for(unsigned int i=0; i<vals.size(); i++){
       map<string,double> &row = vals[i];

       // put the vector of offsets for the current ring into the offsets vector
       if (straw_count==numstraws[ring_count]){
          straw_count=0;
          ring_count++;

          cdc_offsets.push_back(tempvec);

          tempvec.clear();
       }

       // Get the offsets from the calibration database 
       cdc_offset_t temp;
       temp.dx_u=row["dxu"];
       //temp.dx_u=0.;

       temp.dy_u=row["dyu"];
       //temp.dy_u=0.;

       temp.dx_d=row["dxd"];
       //temp.dx_d=0.;

       temp.dy_d=row["dyd"];
       //temp.dy_d=0.;

       tempvec.push_back(temp);

       straw_count++;
    }
    cdc_offsets.push_back(tempvec);
 }
 else{
    jerr<< "CDC wire alignment table not available... bailing... " <<endl;
    exit(0);
 }

 // First axial layer
 for (unsigned int ring=1;ring<5;ring++){ 
    vector<DCDCWire*>straws;
    if (!GetCDCAxialWires(ring,numstraws[ring-1],zcenter,L,cdc_offsets,straws,
             rot_angles,cdc_origin[0],cdc_origin[1])) return false;    
    cdcwires.push_back(straws);
 }  

 // First set of stereo layers
 for (unsigned int i=0;i<8;i++){
    vector<DCDCWire*>straws;
    if (!GetCDCStereoWires(i+5,numstraws[i+4],zcenter,L,cdc_offsets,straws,
             rot_angles,cdc_origin[0],cdc_origin[1])) return false;
    cdcwires.push_back(straws);
 }

 // Second axial layer
 for (unsigned int ring=13;ring<17;ring++){ 
    vector<DCDCWire*>straws;
    if (!GetCDCAxialWires(ring,numstraws[ring-1],zcenter,L,cdc_offsets,straws,
             rot_angles,cdc_origin[0],cdc_origin[1])) return false;    
    cdcwires.push_back(straws);
 }

 // Second set of stereo layers
 for (unsigned int i=8;i<16;i++){
    vector<DCDCWire*>straws;
    if (!GetCDCStereoWires(i+9,numstraws[i+8],zcenter,L,cdc_offsets,straws,
             rot_angles,cdc_origin[0],cdc_origin[1])) return false;
    cdcwires.push_back(straws);
 }

 // Third axial layer
 for (unsigned int ring=25;ring<29;ring++){ 
    vector<DCDCWire*>straws;
    if (!GetCDCAxialWires(ring,numstraws[ring-1],zcenter,L,cdc_offsets,straws,
             rot_angles,cdc_origin[0],cdc_origin[1])) return false;    
    cdcwires.push_back(straws);
 }

 // Calculate wire lengths and compute "s" and "t" direction vectors (orthogonal to "u")
 for (unsigned int i=0;i<cdcwires.size();i++){ 
    for (unsigned int j=0;j<cdcwires[i].size();j++){
       DCDCWire *w=cdcwires[i][j];
       w->L=L/cos(w->stereo);

       // With the addition of close-packed stereo wires, the vector connecting
       // the center of the wire to the beamline ("s" direction) is not necessarily
       // perpendicular to the beamline. By definition, we want the "s" direction
       // to be perpendicular to the wire direction "u" and pointing at the beamline.
       // 
       // NOTE: This extensive comment is here because the result, when implmented
       // below caused a WORSE residual distribution in the close-packed stereo
       // layers. Owing to lack of time currently to track the issue down (most
       // likely in DReferenceTrajectory) I'm commenting out the "correct" calculation
       // of s, but leaving this comment so the issue can be revisited later. This
       // error leads to around 100 micron errors in the C.P.S. wires, but they
       // are completely washed out when the position smearing of 150 microns is applied
       // making the error unnoticable except when position smearing is not applied.
       //
       // April 2, 2009  D.L.
       //
       // Here is how this is calculated -- We define a vector equation with 2 unknowns
       // Z and S:
       //
       //    Zz + Ss = W
       //
       // where:  z = unit vector in z direction
       //         s = unit vector in "s" direction
       //         W = vector pointing to center of wire in lab coordinates
       //
       //  Rearranging, we get:
       //
       //     s = (W - Zz)/S
       //
       //  Since s must be perpendicular to u, we take a dot product of s and u
       // and set it equal to zero to determine Z:
       //
       //    u.s = 0 = u.(W - Zz)/S  =>  u.W = Zu.z
       //
       //   or
       //
       //     Z = u.W/u.z
       //
       //  Thus, the s direction is just given by (W - (u.W/u.z)z)
       //

       //w->sdir=w->origin-DVector3(0,0,w->origin.Z());
       w->sdir = w->origin - DVector3(0.0, 0.0, w->udir.Dot(w->origin)/w->udir.Z());  // see above comments
       w->sdir.SetMag(1.0);

       w->tdir = w->udir.Cross(w->sdir);
       w->tdir.SetMag(1.0); // This isn't really needed
    }
 }  

 return true;
1040}


1043//---------------------------------
1044// GetFDCCathodes
1045//---------------------------------
1046bool DGeometry::GetFDCCathodes(vector<vector<DFDCCathode *> >&fdccathodes) const{
 // Get offsets tweaking nominal geometry from calibration database
 JCalibration * jcalib = dapp->GetJCalibration(runnumber);
 vector<map<string,double> >vals;
 vector<fdc_cathode_offset_t>fdc_cathode_offsets;
 if (jcalib->Get("FDC/cathode_alignment",vals)==false){
    for(unsigned int i=0; i<vals.size(); i++){
       map<string,double> &row = vals[i];

       // Get the offsets from the calibration database 
       fdc_cathode_offset_t temp;
       temp.du=row["dU"];
       //temp.du=0.;

       temp.dphi=row["dPhiU"];
       //temp.dphi=0.;

       fdc_cathode_offsets.push_back(temp);

       temp.du=row["dV"];
       //temp.du=0.;

       temp.dphi=row["dPhiV"];
       //temp.dphi=0.;

       fdc_cathode_offsets.push_back(temp);
    }
 }
 vector< vector<double> >fdc_cathode_pitches;
 if (jcalib->Get("FDC/strip_pitches_v2",vals)==false){
    for(unsigned int i=0; i<vals.size(); i++){
       map<string,double> &row = vals[i];

       vector<double> uvals;
       // Get the offsets from the calibration database 
       uvals.push_back(row["U_SP_1"]);
       uvals.push_back(row["U_G_1"]);
       uvals.push_back(row["U_SP_2"]);
       uvals.push_back(row["U_G_2"]);
       uvals.push_back(row["U_SP_3"]);

       fdc_cathode_pitches.push_back(uvals);

       vector<double> vvals;
       // Get the offsets from the calibration database
       vvals.push_back(row["V_SP_1"]);
       vvals.push_back(row["V_G_1"]);
       vvals.push_back(row["V_SP_2"]);
       vvals.push_back(row["V_G_2"]);
       vvals.push_back(row["V_SP_3"]);

       fdc_cathode_pitches.push_back(vvals);
    }
 }
 else{
    jerr << "Strip pitch calibration unavailable -- setting default..." <<endl;
    // set some sensible default
    for (unsigned int i=0;i<2*FDC_NUM_LAYERS24;i++){

       vector<double> val;
       for (int j = 0; j < 5; j++){
          val.push_back(STRIP_SPACING0.500);
       }

       fdc_cathode_pitches.push_back(val);
    }
 }

 // Generate the vector of cathode plane parameters
 for (int i=0;i<2*FDC_NUM_LAYERS24; i++){
    double angle=(i%2)?(M_PI3.14159265358979323846-CATHODE_ROT_ANGLE1.309):(CATHODE_ROT_ANGLE1.309);

    angle+=fdc_cathode_offsets[i].dphi;
    double SP1 = fdc_cathode_pitches[i][0];
    double SG1 = fdc_cathode_pitches[i][1];
    double SP2 = fdc_cathode_pitches[i][2];
    double SG2 = fdc_cathode_pitches[i][3];
    double SP3 = fdc_cathode_pitches[i][4];

    vector<DFDCCathode *>temp;
    for (int j=0; j<STRIPS_PER_PLANE192; j++){
       DFDCCathode *c = new DFDCCathode;
       c->layer = i+1;
       c->strip = j+1;
       c->angle = angle;
       if (j<48) c->u=(-47.5*SP2 - SG1 + (j-47)*SP1) + fdc_cathode_offsets[i].du;
       else if (j<144) c->u=(double(j)-95.5)*SP2 + fdc_cathode_offsets[i].du;
       else c->u=(47.5*SP2 + SG2 + (j-144)*SP3) + fdc_cathode_offsets[i].du;

       temp.push_back(c);
    }
    fdccathodes.push_back(temp);
 }

 return true;
1141}

1143//---------------------------------
1144// GetFDCWires
1145//---------------------------------
1146bool DGeometry::GetFDCWires(vector<vector<DFDCWire *> >&fdcwires) const{
 // Get geometrical information from database
 vector<double>z_wires;
 vector<double>stereo_angles;

 if(!GetFDCZ(z_wires)) return false;
 if(!GetFDCStereo(stereo_angles)) return false;

 // Get package rotation angles
 double ThetaX[4],ThetaY[4],ThetaZ[4];
 vector<double>rot_angles;
 Get("//posXYZ[@volume='forwardDC_package_1']/@rot", rot_angles);
 ThetaX[0]=rot_angles[0]*M_PI3.14159265358979323846/180.;
 ThetaY[0]=rot_angles[1]*M_PI3.14159265358979323846/180.;
 ThetaZ[0]=rot_angles[2]*M_PI3.14159265358979323846/180.; 
 Get("//posXYZ[@volume='forwardDC_package_2']/@rot", rot_angles);
 ThetaX[1]=rot_angles[0]*M_PI3.14159265358979323846/180.;
 ThetaY[1]=rot_angles[1]*M_PI3.14159265358979323846/180.;
 ThetaZ[1]=rot_angles[2]*M_PI3.14159265358979323846/180.;  
 Get("//posXYZ[@volume='forwardDC_package_3']/@rot", rot_angles);
 ThetaX[2]=rot_angles[0]*M_PI3.14159265358979323846/180.;
 ThetaY[2]=rot_angles[1]*M_PI3.14159265358979323846/180.;
 ThetaZ[2]=rot_angles[2]*M_PI3.14159265358979323846/180.; 
 Get("//posXYZ[@volume='forwardDC_package_4']/@rot", rot_angles);
 ThetaX[3]=rot_angles[0]*M_PI3.14159265358979323846/180.;
 ThetaY[3]=rot_angles[1]*M_PI3.14159265358979323846/180.;
 ThetaZ[3]=rot_angles[2]*M_PI3.14159265358979323846/180.;
 // Get package offsets
 double dX[4],dY[4];
 vector<double>offsets;
 Get("//posXYZ[@volume='forwardDC_package_1']/@X_Y_Z",offsets);
 dX[0]=offsets[0];
 dY[0]=offsets[1];
 Get("//posXYZ[@volume='forwardDC_package_2']/@X_Y_Z",offsets);
 dX[1]=offsets[0];
 dY[1]=offsets[1]; 
 Get("//posXYZ[@volume='forwardDC_package_3']/@X_Y_Z",offsets);
 dX[2]=offsets[0];
 dY[2]=offsets[1];
 Get("//posXYZ[@volume='forwardDC_package_4']/@X_Y_Z",offsets);
 dX[3]=offsets[0];
 dY[3]=offsets[1];

 // Get offsets tweaking nominal geometry from calibration database
 JCalibration * jcalib = dapp->GetJCalibration(runnumber);
 vector<map<string,double> >vals;
 vector<fdc_wire_offset_t>fdc_wire_offsets;
 if (jcalib->Get("FDC/wire_alignment",vals)==false){
    for(unsigned int i=0; i<vals.size(); i++){
       map<string,double> &row = vals[i];

       // Get the offsets from the calibration database 
       fdc_wire_offset_t temp;
       temp.du=row["dU"];
       //temp.du=0.;

       temp.dphi=row["dPhi"];
       //temp.dphi=0.;

       temp.dz=row["dZ"];
       //  temp.dz=0.;

       fdc_wire_offsets.push_back(temp);
    }
 }

 vector<fdc_wire_rotation_t>fdc_wire_rotations;
 if (jcalib->Get("FDC/cell_rotations",vals)==false){
    for(unsigned int i=0; i<vals.size(); i++){
       map<string,double> &row = vals[i];

       // Get the offsets from the calibration database
       fdc_wire_rotation_t temp;
       temp.dPhiX=row["dPhiX"];
       temp.dPhiY=row["dPhiY"];
       temp.dPhiZ=row["dPhiZ"];

       fdc_wire_rotations.push_back(temp);
    }
 }

 // Generate the vector of wire plane parameters
 for(int i=0; i<FDC_NUM_LAYERS24; i++){
    double angle=-stereo_angles[i]*M_PI3.14159265358979323846/180.+fdc_wire_offsets[i].dphi;

    vector<DFDCWire *>temp;
    for(int j=0; j<WIRES_PER_PLANE96; j++){
       unsigned int pack_id=i/6;

       DFDCWire *w = new DFDCWire;
       w->layer = i+1;
       w->wire = j+1;
       w->angle = angle;

       // find coordinates of center of wire in rotated system
       float u = U_OF_WIRE_ZERO(-((96 -1)*1.0)/2) + WIRE_SPACING1.0*(float)(j);
       w->u=u+fdc_wire_offsets[i].du;

       // Rotate coordinates into lab system and set the wire's origin
       // Note that the FDC measures "angle" such that angle=0
       // corresponds to the anode wire in the vertical direction
       // (i.e. at phi=90 degrees).
       float x = u*sin(angle + M_PI3.14159265358979323846/2.0);
       float y = u*cos(angle + M_PI3.14159265358979323846/2.0);
       w->origin.SetXYZ(x,y,0.);
       w->origin.RotateX(ThetaX[pack_id]+fdc_wire_rotations[i].dPhiX);
       w->origin.RotateY(ThetaY[pack_id]+fdc_wire_rotations[i].dPhiY);
       w->origin.RotateZ(ThetaZ[pack_id]+fdc_wire_rotations[i].dPhiZ);
       DVector3 globalOffsets(dX[pack_id],dY[pack_id],z_wires[i]+fdc_wire_offsets[i].dz);
       w->origin+=globalOffsets;

       // Length of wire is set by active radius
       w->L = 2.0*sqrt(pow(FDC_ACTIVE_RADIUS48.5,2.0) - u*u);

       // Set directions of wire's coordinate system with "udir"
       // along wire.
       w->udir.SetXYZ(sin(angle),cos(angle),0.0);
       w->udir.RotateX(ThetaX[pack_id]+fdc_wire_rotations[i].dPhiX);
       w->udir.RotateY(ThetaY[pack_id]+fdc_wire_rotations[i].dPhiY);
       w->udir.RotateZ(ThetaZ[pack_id]+fdc_wire_rotations[i].dPhiZ);
       w->angles.SetXYZ(ThetaX[pack_id]+fdc_wire_rotations[i].dPhiX,
             ThetaY[pack_id]+fdc_wire_rotations[i].dPhiY,
             ThetaZ[pack_id]+fdc_wire_rotations[i].dPhiZ);
       w->u+=dX[pack_id]*w->udir.y()-dY[pack_id]*w->udir.x();

       // "s" points in direction from beamline to midpoint of
       // wire. This happens to be the same direction as "origin"
       w->sdir = w->origin;
       w->sdir.SetMag(1.0);

       w->tdir = w->udir.Cross(w->sdir);
       w->tdir.SetMag(1.0); // This isn't really needed
       temp.push_back(w);
    }
    fdcwires.push_back(temp);
 }

 return true;
1284}

1286//---------------------------------
1287// GetFDCZ
1288//---------------------------------
1289bool DGeometry::GetFDCZ(vector<double> &z_wires) const
1290{
 // The FDC geometry is defined as 4 packages, each containing 2
 // "module"s and each of those containing 3 "chambers". The modules
 // are placed as multiple copies in Z using mposZ, but none of the
 // others are (???).
 //
 // This method is currently hardwired to assume 4 packages and
 // 3 chambers. (The number of modules is discovered via the
 // "ncopy" attribute of mposZ.)

 vector<double> ForwardDC;
 vector<double> forwardDC;
 vector<double> forwardDC_package[4];
 vector<double> forwardDC_module[4];
 vector<double> forwardDC_chamber[4][6];

 if(!Get("//section/composition/posXYZ[@volume='ForwardDC']/@X_Y_Z", ForwardDC)) return false;
 if(!Get("//composition[@name='ForwardDC']/posXYZ[@volume='forwardDC']/@X_Y_Z", forwardDC)) return false;
 if(!Get("//posXYZ[@volume='forwardDC_package_1']/@X_Y_Z", forwardDC_package[0])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_package_2']/@X_Y_Z", forwardDC_package[1])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_package_3']/@X_Y_Z", forwardDC_package[2])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_package_4']/@X_Y_Z", forwardDC_package[3])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_module_1']/@X_Y_Z", forwardDC_module[0])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_module_2']/@X_Y_Z", forwardDC_module[1])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_module_3']/@X_Y_Z", forwardDC_module[2])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_module_4']/@X_Y_Z", forwardDC_module[3])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='1']", forwardDC_chamber[0][0])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='2']", forwardDC_chamber[0][1])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='3']", forwardDC_chamber[0][2])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='4']", forwardDC_chamber[0][3])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='5']", forwardDC_chamber[0][4])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='6']", forwardDC_chamber[0][5])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='1']", forwardDC_chamber[1][0])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='2']", forwardDC_chamber[1][1])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='3']", forwardDC_chamber[1][2])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='4']", forwardDC_chamber[1][3])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='5']", forwardDC_chamber[1][4])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='6']", forwardDC_chamber[1][5])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='1']", forwardDC_chamber[2][0])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='2']", forwardDC_chamber[2][1])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='3']", forwardDC_chamber[2][2])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='4']", forwardDC_chamber[2][3])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='5']", forwardDC_chamber[2][4])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='6']", forwardDC_chamber[2][5])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='1']", forwardDC_chamber[3][0])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='2']", forwardDC_chamber[3][1])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='3']", forwardDC_chamber[3][2])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='4']", forwardDC_chamber[3][3])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='5']", forwardDC_chamber[3][4])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='6']", forwardDC_chamber[3][5])) return false;

 // Offset due to global FDC envelopes
 double zfdc = ForwardDC[2] + forwardDC[2];

 // Loop over packages
 for(int package=1; package<=4; package++){
    double z_package = forwardDC_package[package-1][2];

    // Each "package" has 1 "module" which is currently positioned at 0,0,0 but
    // that could in principle change so we add the module z-offset
    double z_module = forwardDC_module[package-1][2];

    // Loop over chambers in this module
    for(int chamber=1; chamber<=6; chamber++){
       double z_chamber = forwardDC_chamber[package-1][chamber-1][2];

       double z = zfdc + z_package + z_module + z_chamber;				
       z_wires.push_back(z);
    }
 }

 return true;
1362}

1364//---------------------------------
1365// GetFDCStereo
1366//---------------------------------
1367bool DGeometry::GetFDCStereo(vector<double> &stereo_angles) const
1368{
 // The FDC geometry is defined as 4 packages, each containing 2
 // "module"s and each of those containing 3 "chambers". The modules
 // are placed as multiple copies in Z using mposZ, but none of the
 // others are (???).
 //
 // This method is currently hardwired to assume 4 packages and
 // 3 chambers. (The number of modules is discovered via the
 // "ncopy" attribute of mposZ.)
 //
 // Stereo angles are assumed to be rotated purely about the z-axis
 // and the units are not specified, but the XML currently uses degrees.

 vector<double> forwardDC_module[4];
 vector<double> forwardDC_chamber[4][6];

 if(!Get("//posXYZ[@volume='forwardDC_module_1']/@X_Y_Z", forwardDC_module[0])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_module_2']/@X_Y_Z", forwardDC_module[1])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_module_3']/@X_Y_Z", forwardDC_module[2])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_module_4']/@X_Y_Z", forwardDC_module[3])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='1']", forwardDC_chamber[0][0])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='2']", forwardDC_chamber[0][1])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='3']", forwardDC_chamber[0][2])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='4']", forwardDC_chamber[0][3])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='5']", forwardDC_chamber[0][4])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='6']", forwardDC_chamber[0][5])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='1']", forwardDC_chamber[1][0])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='2']", forwardDC_chamber[1][1])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='3']", forwardDC_chamber[1][2])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='4']", forwardDC_chamber[1][3])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='5']", forwardDC_chamber[1][4])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='6']", forwardDC_chamber[1][5])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='1']", forwardDC_chamber[2][0])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='2']", forwardDC_chamber[2][1])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='3']", forwardDC_chamber[2][2])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='4']", forwardDC_chamber[2][3])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='5']", forwardDC_chamber[2][4])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='6']", forwardDC_chamber[2][5])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='1']", forwardDC_chamber[3][0])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='2']", forwardDC_chamber[3][1])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='3']", forwardDC_chamber[3][2])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='4']", forwardDC_chamber[3][3])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='5']", forwardDC_chamber[3][4])) return false;
 if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='6']", forwardDC_chamber[3][5])) return false;

 // Loop over packages
 for(int package=1; package<=4; package++){

    // Loop over chambers
    for(int chamber=1; chamber<=6; chamber++){
       // if (chamber==4) forwardDC_chamber[package-1][chamber-1][2]+=15.0;
       stereo_angles.push_back(forwardDC_chamber[package-1][chamber-1][2]);
    }
 }

 return true;
1424}

1426//---------------------------------
1427// GetFDCRmin
1428//---------------------------------
1429bool DGeometry::GetFDCRmin(vector<double> &rmin_packages) const
1430{
 vector<double> FDA[4];

 if(!Get("//section[@name='ForwardDC']/tubs[@name='FDA1']/@Rio_Z", FDA[0])) return false;
 if(!Get("//section[@name='ForwardDC']/tubs[@name='FDA2']/@Rio_Z", FDA[1])) return false;
 if(!Get("//section[@name='ForwardDC']/tubs[@name='FDA3']/@Rio_Z", FDA[2])) return false;
 if(!Get("//section[@name='ForwardDC']/tubs[@name='FDA4']/@Rio_Z", FDA[3])) return false;

 rmin_packages.push_back(FDA[0][0]);
 rmin_packages.push_back(FDA[1][0]);
 rmin_packages.push_back(FDA[2][0]);
 rmin_packages.push_back(FDA[3][0]);

 return true;
1444}

1446//---------------------------------
1447// GetFDCRmax
1448//---------------------------------
1449bool DGeometry::GetFDCRmax(double &rmax_active_fdc) const
1450{
 // We assume that all packages have the same outer radius of the
 // active area.
 vector<double> FDA1;

 bool good = Get("//section[@name='ForwardDC']/tubs[@name='FDA1']/@Rio_Z", FDA1);

 if(!good){
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1458<<" "<<"Unable to retrieve FDC Rmax values."<<endl;
    return good;
 }

 rmax_active_fdc = FDA1[1];

 return good;
1465}

1467//---------------------------------
1468// GetCDCOption
1469//---------------------------------
1470bool DGeometry::GetCDCOption(string &cdc_option) const
1471{
 bool good = Get("//section[@name='CentralDC_s']/composition/posXYZ/@volume", cdc_option);

 if(!good){
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1475<<" "<<"Unable to retrieve CDC option string."<<endl;
 }

 return good;
1479}

1481//---------------------------------
1482// GetCDCCenterZ
1483//---------------------------------
1484bool DGeometry::GetCDCCenterZ(double &cdc_center_z) const
1485{

 return false;
1488}

1490//---------------------------------
1491// GetCDCAxialLength
1492//---------------------------------
1493bool DGeometry::GetCDCAxialLength(double &cdc_axial_length) const
1494{
 vector<double> Rio_Z;
 bool good = Get("//section[@name='CentralDC']/tubs[@name='STRW']/@Rio_Z", Rio_Z);
 cdc_axial_length = Rio_Z[2];

 if(!good){
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1500<<" "<<"Unable to retrieve CDC axial wire length"<<endl;
 }

 return false;
1504}

1506//---------------------------------
1507// GetCDCStereo
1508//---------------------------------
1509bool DGeometry::GetCDCStereo(vector<double> &cdc_stereo) const
1510{

 return false;
1513}

1515//---------------------------------
1516// GetCDCRmid
1517//---------------------------------
1518bool DGeometry::GetCDCRmid(vector<double> &cdc_rmid) const
1519{

 return false;
1522}

1524//---------------------------------
1525// GetCDCNwires
1526//---------------------------------
1527bool DGeometry::GetCDCNwires(vector<int> &cdc_nwires) const
1528{
cdc_nwires.push_back(42);
cdc_nwires.push_back(42);
cdc_nwires.push_back(54);
cdc_nwires.push_back(54);
cdc_nwires.push_back(66);
cdc_nwires.push_back(66);
cdc_nwires.push_back(80);
cdc_nwires.push_back(80);
cdc_nwires.push_back(93);
cdc_nwires.push_back(93);
cdc_nwires.push_back(106);
cdc_nwires.push_back(106);
cdc_nwires.push_back(123);
cdc_nwires.push_back(123);
cdc_nwires.push_back(135);
cdc_nwires.push_back(135);
cdc_nwires.push_back(146);
cdc_nwires.push_back(146);
cdc_nwires.push_back(158);
cdc_nwires.push_back(158);
cdc_nwires.push_back(170);
cdc_nwires.push_back(170);
cdc_nwires.push_back(182);
cdc_nwires.push_back(182);
cdc_nwires.push_back(197);
cdc_nwires.push_back(197);
cdc_nwires.push_back(209);
cdc_nwires.push_back(209);

 return false;
1559}


1562//---------------------------------
1563// GetCDCEndplate
1564//---------------------------------
1565bool DGeometry::GetCDCEndplate(double &z,double &dz,double &rmin,double &rmax)
 const{

    vector<double>cdc_origin;
    vector<double>cdc_center;
    vector<double>cdc_layers_offset;
    vector<double>cdc_endplate_pos;
    vector<double>cdc_endplate_dim;

    if(!Get("//posXYZ[@volume='CentralDC'/@X_Y_Z",cdc_origin)) return false;
    if(!Get("//posXYZ[@volume='centralDC']/@X_Y_Z",cdc_center)) return false;
    if(!Get("//posXYZ[@volume='CDPD']/@X_Y_Z",cdc_endplate_pos)) return false;
    if(!Get("//tubs[@name='CDPD']/@Rio_Z",cdc_endplate_dim)) return false;
    if(!Get("//posXYZ[@volume='CDClayers']/@X_Y_Z",cdc_layers_offset)) return false;

    if(cdc_origin.size()<3){
       _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1581<<" "<<"cdc_origin.size()<3 !"<<endl;
       return false;
    }
    if(cdc_center.size()<3){
       _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1585<<" "<<"cdc_center.size()<3 !"<<endl;
       return false;
    }
    if(cdc_endplate_pos.size()<3){
       _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1589<<" "<<"cdc_endplate_pos.size()<3 !"<<endl;
       return false;
    }
    if(cdc_endplate_dim.size()<3){
       _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1593<<" "<<"cdc_endplate_dim.size()<3 !"<<endl;
       return false;
    }
    if (cdc_layers_offset.size()<3){
       _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1597<<" "<<"cdc_layers_offset.size()<3 !"<<endl;
       return false;
    }

    z=cdc_origin[2]+cdc_center[2]+cdc_endplate_pos[2]+cdc_layers_offset[2];
    dz=cdc_endplate_dim[2];
    rmin=cdc_endplate_dim[0];
    rmax=cdc_endplate_dim[1];

    return true;
 }
1608//---------------------------------
1609// GetBCALRmin
1610//---------------------------------
1611// Including the support plate
1612bool DGeometry::GetBCALRmin(float &bcal_rmin) const
1613{
 vector<float> bcal_mother_Rio_Z;
 bool good = Get("//section[@name='BarrelEMcal']/tubs[@name='BCAL']/@Rio_Z", bcal_mother_Rio_Z);
 if(!good){
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1617<<" "<<"Unable to retrieve BCAL mother RioZ info."<<endl;
    bcal_rmin = 0.0;
    return false;
 }
 if(bcal_mother_Rio_Z.size() == 3){
    bcal_rmin = bcal_mother_Rio_Z[0];
    return true;
 }
 else{
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1626<<" "<<"Wrong vector size for BCAL mother RioZ!!!"<<endl;
    bcal_rmin = 0.0;
    return false;
 }
1630}

1632//---------------------------------
1633// GetBCALfADCRadii
1634//---------------------------------
1635bool DGeometry::GetBCALfADCRadii(vector<float> &fADC_radii) const
1636{
 vector<float> BM[5];

 if(!Get("//section[@name='BarrelEMcal']/tubs[@name='BM01']/@Rio_Z", BM[0])) return false;
 if(!Get("//section[@name='BarrelEMcal']/tubs[@name='BM02']/@Rio_Z", BM[1])) return false;
 if(!Get("//section[@name='BarrelEMcal']/tubs[@name='BM04']/@Rio_Z", BM[2])) return false;
 if(!Get("//section[@name='BarrelEMcal']/tubs[@name='BMF7']/@Rio_Z", BM[3])) return false;
 if(!Get("//section[@name='BarrelEMcal']/tubs[@name='BMFA']/@Rio_Z", BM[4])) return false;

 fADC_radii.push_back(BM[0][0]);
 fADC_radii.push_back(BM[1][0]);
 fADC_radii.push_back(BM[2][0]);
 fADC_radii.push_back(BM[3][0]);
 fADC_radii.push_back(BM[4][1]);

 return true;
1652}

1654//---------------------------------
1655// GetBCALNmodules
1656//---------------------------------
1657bool DGeometry::GetBCALNmodules(unsigned int &bcal_nmodules) const
1658{
 vector<unsigned int> ncopy;
 bool good = Get("//section[@name='BarrelEMcal']/composition/mposPhi/@ncopy", ncopy);
 if(!good){
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1662<<" "<<"Unable to retrieve BCAL barrelModule ncopy info."<<endl;
    bcal_nmodules = 0;
    return false;
 }
 if(ncopy.size() == 1){
    bcal_nmodules = ncopy[0];
    return true;
 }else{
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1670<<" "<<"Wrong vector size for BCAL barrelModule ncopy!!!"<<endl;
    bcal_nmodules = 0;
    return false;
 }
1674}

1676//---------------------------------
1677// GetBCALCenterZ
1678//---------------------------------
1679bool DGeometry::GetBCALCenterZ(float &bcal_center_z) const
1680{
 vector<float> z0;
 bool good = Get("//section[@name='BarrelEMcal']/parameters/real[@name='z0']/@value", z0);
 if(!good){
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1684<<" "<<"Unable to retrieve BCAL parameters z0 info."<<endl;
    bcal_center_z = 0.0;
    return false;
 }
 if(z0.size() == 1){
    bcal_center_z = z0[0];
    return true;
 }else{
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1692<<" "<<"Wrong vector size for BCAL parameters z0!!!"<<endl;
    bcal_center_z = 0.0;
    return false;
 }

1697}

1699//---------------------------------
1700// GetBCALLength
1701//---------------------------------
1702// The lightguides are not included
1703bool DGeometry::GetBCALLength(float &bcal_length) const
1704{
 vector<float> module_length;
 bool good = Get("//section[@name='BarrelEMcal']/tubs[@name='BM01']/@Rio_Z", module_length);
 if(!good){
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1708<<" "<<"Unable to retrieve BCAL submodule RioZ info."<<endl;
    bcal_length = 0.0;
    return false;
 }
 if(module_length.size() == 3){
    bcal_length = module_length[2];
    return true;
 }
 else{
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1717<<" "<<"Wrong vector size for BCAL submodule RioZ!!!"<<endl;
    bcal_length = 0.0;
    return false;
 }
1721}

1723//---------------------------------
1724// GetBCALDepth
1725//---------------------------------
1726// Including the support plate and the support bar
1727bool DGeometry::GetBCALDepth(float &bcal_depth) const
1728{
 vector<float> bcal_moth_Rio_Z;
 bool good = Get("//section[@name='BarrelEMcal']/tubs[@name='BCAL']/@Rio_Z", bcal_moth_Rio_Z);
 if(!good){
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1732<<" "<<"Unable to retrieve BCAL mother RioZ info."<<endl;
    bcal_depth = 0.0;
    return false;
 }
 if(bcal_moth_Rio_Z.size() == 3){
    bcal_depth = bcal_moth_Rio_Z[1] - bcal_moth_Rio_Z[0];
    return true;
 }
 else{
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1741<<" "<<"Wrong vector size for BCAL mother RioZ!!!"<<endl;
    bcal_depth = 0.0;
    return false;
 }

1746}

1748//---------------------------------
1749// GetBCALPhiShift
1750//---------------------------------
1751bool DGeometry::GetBCALPhiShift(float &bcal_phi_shift) const
1752{
 vector<float> Phi0;
 bool good = Get("//section[@name='BarrelEMcal']/composition/mposPhi/@Phi0", Phi0);
 if(!good) return false;
 if(Phi0.size() == 1){
    bcal_phi_shift = Phi0[0];
    return true;
 }else{
    bcal_phi_shift = 0.0;
    return false;
 }
1763}

1765//---------------------------------
1766// GetCCALZ
1767//---------------------------------
1768bool DGeometry::GetCCALZ(double &z_ccal) const
1769{
 vector<double> ComptonEMcalpos;
 jgeom->SetVerbose(0);   // don't print error messages for optional detector elements
 bool good = Get("//section/composition/posXYZ[@volume='ComptonEMcal']/@X_Y_Z", ComptonEMcalpos);\
 jgeom->SetVerbose(1);   // reenable error messages

 if(!good){
 // NEED TO RETHINK ERROR REPORTING FOR OPTIONAL DETECTOR ELEMENTS
    //_DBG_<<"Unable to retrieve ComptonEMcal position."<<endl;  
    z_ccal = 1279.376;   
    return false;
 }else{
  z_ccal = ComptonEMcalpos[2];
    return true;
 }
1784}

1786//---------------------------------
1787// GetFCALZ
1788//---------------------------------
1789bool DGeometry::GetFCALZ(double &z_fcal) const
1790{
 vector<double> ForwardEMcalpos;
 bool good = Get("//section/composition/posXYZ[@volume='ForwardEMcal']/@X_Y_Z", ForwardEMcalpos);

 if(!good){
    _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1795<<" "<<"Unable to retrieve ForwardEMcal position."<<endl;
    z_fcal=0.0;
    return false;
 }else{
    z_fcal = ForwardEMcalpos[2];
    return true;
 }
1802}


1805//---------------------------------
1806// GetFCALPosition
1807//---------------------------------
1808bool DGeometry::GetFCALPosition(double &x,double &y,double &z) const
1809{
vector<double> ForwardEMcalpos;
bool good = Get("//section/composition/posXYZ[@volume='ForwardEMcal']/@X_Y_Z", ForwardEMcalpos);

if(!good){
  _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":"
<<1814<<" "<<"Unable to retrieve ForwardEMcal position."<<endl;
  x=0.,y=0.,z=0.;
  return false;
}else{
  x=ForwardEMcalpos[0],y=ForwardEMcalpos[1],z=ForwardEMcalpos[2];
  //_DBG_ << "FCAL position: (x,y,z)=(" << x <<"," << y << "," << z << ")"<<endl;
  return true;
}
1822}

1824//---------------------------------
1825// GetCCALPosition
1826//---------------------------------
1827bool DGeometry::GetCCALPosition(double &x,double &y,double &z) const
1828{
vector<double> ComptonEMcalpos;
jgeom->SetVerbose(0);   // don't print error messages for optional detector elements
bool good = Get("//section/composition/posXYZ[@volume='ComptonEMcal']/@X_Y_Z", ComptonEMcalpos);
jgeom->SetVerbose(1);   // reenable error messages

if(!good){
  //_DBG_<<"Unable to retrieve ComptonEMcal position."<<endl;
  x=0.,y=0.,z=0.;
  return false;
}else{
  x=ComptonEMcalpos[0],y=ComptonEMcalpos[1],z=ComptonEMcalpos[2]; 
  //_DBG_ << "CCAL position: (x,y,z)=(" << ComptonEMcalpos[0] <<","
// << ComptonEMcalpos[1]<<","<<ComptonEMcalpos[2]<< ")" << endl;
  return true;
}
1844}

1846//---------------------------------
1847// GetFCALInsertRowSize
1848//---------------------------------
1849bool DGeometry::GetFCALInsertRowSize(int &insert_row_size) const
1850{
 jgeom->SetVerbose(0);   // don't print error messages for optional detector elements
 bool good = Get("//composition[@name='LeadTungstateFullRow']/mposX[@volume='LTBLwrapped']/@ncopy",insert_row_size);
 jgeom->SetVerbose(1);   // reenable error messages

 if(!good){
 // NEED TO RETHINK ERROR REPORTING FOR OPTIONAL DETECTOR ELEMENTS
    //_DBG_<<"Unable to retrieve ComptonEMcal position."<<endl;  
    insert_row_size = 0;   
    return false;
 }else{
    return true;
 }
1863}

1865//---------------------------------
1866// GetFCALBlockSize
1867//---------------------------------
1868bool DGeometry::GetFCALBlockSize(vector<double> &block) const
1869{
 jgeom->SetVerbose(0);   // don't print error messages for optional detector elements
 bool good = Get("//box[@name='LGBL']/@X_Y_Z",block);
 jgeom->SetVerbose(1);   // reenable error messages

 if(!good){
 // NEED TO RETHINK ERROR REPORTING FOR OPTIONAL DETECTOR ELEMENTS
    //_DBG_<<"Unable to retrieve ComptonEMcal position."<<endl;  
    return false;
 }else{
    return true;
 }
1881}

1883//---------------------------------
1884// GetFCALInsertBlockSize
1885//---------------------------------
1886bool DGeometry::GetFCALInsertBlockSize(vector<double> &block) const
1887{
 jgeom->SetVerbose(0);   // don't print error messages for optional detector elements
 bool good = Get("//box[@name='LTB1']/@X_Y_Z",block);
 jgeom->SetVerbose(1);   // reenable error messages

 if(!good){
 // NEED TO RETHINK ERROR REPORTING FOR OPTIONAL DETECTOR ELEMENTS
    //_DBG_<<"Unable to retrieve ComptonEMcal position."<<endl;  
    return false;
 }else{
    return true;
 }
1899}

1901//---------------------------------
1902// GetDIRCZ
1903//---------------------------------
1904bool DGeometry::GetDIRCZ(double &z_dirc) const
1905{
vector<double> dirc_face;
jgeom->SetVerbose(0);   // don't print error messages for optional detector elements
bool good = Get("//section/composition/posXYZ[@volume='DIRC']/@X_Y_Z",dirc_face);
jgeom->SetVerbose(1);   // reenable error messages

if(!good){
  //_DBG_<<"Unable to retrieve DIRC position."<<endl;
  z_dirc=0.0;
  return false;
}
else{ 
  vector<double>dirc_plane;
  vector<double>dirc_shift;
  vector<double>bar_plane;
  Get("//composition[@name='DRCC']/posXYZ[@volume='DCML10']/@X_Y_Z/plane[@value='1']", dirc_plane);
  Get("//composition[@name='DIRC']/posXYZ[@volume='DRCC']/@X_Y_Z", dirc_shift);
  z_dirc=dirc_face[2]+dirc_plane[2]+dirc_shift[2] + 0.8625; // last shift is the average center of quartz bar (585.862)

  //jout << "DIRC z position = " << z_dirc << " cm." << endl;
  return true;
}
1927}

1929//---------------------------------
1930// GetTRDZ
1931//---------------------------------
1932bool DGeometry::GetTRDZ(vector<double> &z_trd) const
1933{
 vector<double> trd_origin;
 jgeom->SetVerbose(0);   // don't print error messages for optional detector elements
 bool good = Get("//section/composition/posXYZ[@volume='TRDGEM']/@X_Y_Z",trd_origin);
 jgeom->SetVerbose(1);   // reenable error messages

 if(!good){
   //_DBG_<<"Unable to retrieve TRD position."<<endl;
   return false;
 }
 else{ 
   vector<double>trd_G10;
   Get("//composition[@name='TRDGEM']/posXYZ[@volume='TGPK']/@X_Y_Z",trd_G10);

   jout << "TRD z positions = ";
   for(int i=0; i<5; i++) {
     vector<double>trd_plane;
     Get(Form("//composition[@name='TGPK']/posXYZ[@volume='TRDG']/@X_Y_Z/plane[@value='%d']",i),trd_plane);
     double z_trd_plane=trd_origin[2]+trd_G10[2]+trd_plane[2];
     jout << z_trd_plane << ", ";
     z_trd.push_back(z_trd_plane);
   }
   jout << "cm" << endl;
 }

 return true;
1959}

1961//---------------------------------
1962// GetTOFZ
1963//---------------------------------
1964bool DGeometry::GetTOFZ(vector<double> &z_tof) const
1965{
 vector<double> ForwardTOF;
 vector<double> forwardTOF[2];
 vector<double> FTOC;

 if(!Get("//section/composition/posXYZ[@volume='ForwardTOF']/@X_Y_Z", ForwardTOF)) return false;
 if(!Get("//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='0']", forwardTOF[0])) return false;
 if(!Get("//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='1']", forwardTOF[1])) return false;
 if(!Get("//box[@name='FTOC' and sensitive='true']/@X_Y_Z", FTOC)) return false;

 z_tof.push_back(ForwardTOF[2] + forwardTOF[0][2] - FTOC[2]/2.0);
 z_tof.push_back(ForwardTOF[2] + forwardTOF[1][2] - FTOC[2]/2.0);

//cerr << "DGeometry::GetTOFZ() = " << z_tof[0] << " " << z_tof[1] << endl;

 return true;
1981}

1983//---------------------------------
1984// GetTOFPaddleParameters
1985//---------------------------------
1986bool DGeometry::GetTOFPaddleParameters(map<string,double> &paddle_params) const
1987{
vector<double> xyz_bar;

// load the number of bars in each area
  int num_bars1 = 0;
  if(!Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOC']/@ncopy",num_bars1)) return false; 
  int num_narrow_bars1 = 0;
  if(!Get("//composition[@name='forwardTOF_bottom2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars1)) return false; 
  int num_single_end_bars1 = 0;
  if(!Get("//composition[@name='forwardTOF_north']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars1)) return false;

  jgeom->SetVerbose(0);   // don't print error messages for optional detector elements
  int num_narrower_bars1 = 0;   // optional - added during upgrade
  Get("//composition[@name='forwardTOF_bottom3']/mposY[@volume='FTOL']/@ncopy",num_narrower_bars1); 
  int num_narrower_bars2 = 0;   // optional  - added during upgrade
  Get("//composition[@name='forwardTOF_top3']/mposY[@volume='FTOL']/@ncopy",num_narrower_bars2); 
  jgeom->SetVerbose(1);   // reenable error messages

  int num_narrow_bars2 = 0;
  if(!Get("//composition[@name='forwardTOF_top2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars2)) return false;
  int num_bars2 = 0;
  if(!Get("//composition[@name='forwardTOF_top1']/mposY[@volume='FTOC']/@ncopy",num_bars2)) return false;
  int num_single_end_bars2 = 0;
  if(!Get("//composition[@name='forwardTOF_south']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars2)) return false;

int NLONGBARS = num_bars1 + num_bars2 + num_narrow_bars1 + num_narrow_bars2
				 + num_narrower_bars1 + num_narrower_bars2;
int NSHORTBARS = num_single_end_bars1 + num_single_end_bars2;
int FIRSTSHORTBAR = num_bars1 + num_narrow_bars1 + num_narrower_bars1 + 1;
int LASTSHORTBAR = FIRSTSHORTBAR + NSHORTBARS/2 - 1;

// load bar sizes
//Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOC']/@X_Y_Z",xyz_bar);
  if(!Get("//box[@name='FTOC' and sensitive='true']/@X_Y_Z", xyz_bar)) return false;
double LONGBARLENGTH = xyz_bar[0];
double BARWIDTH  = xyz_bar[1];
//Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOH']/@X_Y_Z",xyz_bar);
  if(!Get("//box[@name='FTOH' and sensitive='true']/@X_Y_Z", xyz_bar)) return false;
double SHORTBARLENGTH = xyz_bar[0];


// load up the structure containing the parameters for the calling function
paddle_params["NLONGBARS"] = NLONGBARS;
paddle_params["NSHORTBARS"] = NSHORTBARS;
paddle_params["BARWIDTH"] = BARWIDTH;

paddle_params["LONGBARLENGTH"] = LONGBARLENGTH;
paddle_params["HALFLONGBARLENGTH"] = LONGBARLENGTH/2.;
paddle_params["SHORTBARLENGTH"] = SHORTBARLENGTH;
paddle_params["HALFSHORTBARLENGTH"] = SHORTBARLENGTH/2.;

paddle_params["FIRSTSHORTBAR"] = FIRSTSHORTBAR;
paddle_params["LASTSHORTBAR"] = LASTSHORTBAR;

//cout << "In DGeometry::GetTOFPaddleParameters() ..." << endl;
//for(auto el : paddle_params) {
 	//	std::cout << el.first << " " << el.second << endl;
//}	

return true;
2047}


2050//---------------------------------
2051// GetTOFPaddlePerpPositions
2052//---------------------------------
2053bool DGeometry::GetTOFPaddlePerpPositions(vector<double> &y_tof, vector<double> &y_widths) const
2054{
// add in a dummy entry, since we are indexing by paddle number, which starts at 1
// maybe change this some day?
y_tof.push_back(0);
y_widths.push_back(0);

	// Next fill array of bar positions within a plane
	// y_tof[barnumber] gives y position in the center of the bar. [currently barnumber = 1 - 46]
	double y0,dy;
1
'y0' declared without an initial value→

// load the number of bars
  int num_bars=1;   // start counting at 1
  int num_bars1 = 0;
  Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOC']/@ncopy",num_bars1); 
  int num_narrow_bars1 = 0;
  Get("//composition[@name='forwardTOF_bottom2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars1); 
  int num_single_end_bars1 = 0;
  Get("//composition[@name='forwardTOF_north']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars1); 

  jgeom->SetVerbose(0);   // don't print error messages for optional detector elements
  int num_narrower_bars1 = 0;
  Get("//composition[@name='forwardTOF_bottom3']/mposY[@volume='FTOL']/@ncopy",num_narrower_bars1); 
  int num_narrower_bars2 = 0;
  Get("//composition[@name='forwardTOF_top3']/mposY[@volume='FTOL']/@ncopy",num_narrower_bars2); 
  jgeom->SetVerbose(1);   // reenable error messages

  int num_narrow_bars2 = 0;
  Get("//composition[@name='forwardTOF_top2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars2); 
  int num_bars2 = 0;
  Get("//composition[@name='forwardTOF_top1']/mposY[@volume='FTOC']/@ncopy",num_bars2); 
  int num_single_end_bars2 = 0;
  Get("//composition[@name='forwardTOF_south']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars2); 

	// First 19 long bars
	Get("//composition[@name='forwardTOF_bottom1']/mposY/@Y0",y0);
2
←
Calling 'DGeometry::Get'→
9
←
Returning from 'DGeometry::Get'→
	Get("//composition[@name='forwardTOF_bottom1']/mposY/@dY",dy);
	vector<double>tof_bottom1;
	Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_bottom1']/@X_Y_Z",tof_bottom1);  
	for (int k=num_bars;k<num_bars+num_bars1;k++){
10
←
Assuming the condition is true→
11
←
Loop condition is true.  Entering loop body→
  	y_tof.push_back(y0+tof_bottom1[1]+dy*double(k-num_bars));
12
←
The left operand of '+' is a garbage value
  	y_widths.push_back(dy);
	}
	num_bars+=num_bars1;

	// two narrow long bars
	Get("//composition[@name='forwardTOF_bottom2']/mposY/@Y0",y0);
	Get("//composition[@name='forwardTOF_bottom2']/mposY/@dY",dy);
	vector<double>tof_bottom2;
	Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_bottom2']/@X_Y_Z",tof_bottom2);  
	for (int k=num_bars;k<num_bars+num_narrow_bars1;k++){
  	y_tof.push_back(y0+tof_bottom2[1]+dy*double(k-num_bars));
  	y_widths.push_back(dy);
	}
	num_bars+=num_narrow_bars1;
	
	// two narrower long bars - added by upgrade
	if(num_narrower_bars1 > 0) {
Get("//composition[@name='forwardTOF_bottom3']/mposY/@Y0",y0);
Get("//composition[@name='forwardTOF_bottom3']/mposY/@dY",dy);
vector<double>tof_bottom3;
Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_bottom3']/@X_Y_Z",tof_bottom3);  
for (int k=num_bars;k<num_bars+num_narrower_bars1;k++){
	y_tof.push_back(y0+tof_bottom3[1]+dy*double(k-num_bars));
   		y_widths.push_back(dy);
}
num_bars+=num_narrow_bars1;
	}

	// two short wide bars  (4 wide in upgrade)
	Get("//composition[@name='forwardTOF_north']/mposY/@Y0",y0);
	Get("//composition[@name='forwardTOF_north']/mposY/@dY",dy);
	vector<double>tof_north;
	Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_north']/@X_Y_Z",tof_north);  
	for (int k=num_bars;k<num_bars+num_single_end_bars1;k++){
  	y_tof.push_back(y0+tof_north[1]+dy*double(k-num_bars));
  	y_widths.push_back(dy);
	}
	num_bars+=num_single_end_bars1;

	// two narrower long bars - added by upgrade
	if(num_narrower_bars2 > 0) {
Get("//composition[@name='forwardTOF_top3']/mposY/@Y0",y0);
Get("//composition[@name='forwardTOF_top3']/mposY/@dY",dy);
vector<double>tof_top3;
Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_top3']/@X_Y_Z",tof_top3);  
for (int k=num_bars;k<num_bars+num_narrower_bars2;k++){
	y_tof.push_back(y0+tof_top3[1]+dy*double(k-num_bars));
  		y_widths.push_back(dy);
}
num_bars+=num_narrow_bars2;
	}

	// two narrow long bars
	Get("//composition[@name='forwardTOF_top2']/mposY/@Y0",y0);
	Get("//composition[@name='forwardTOF_top2']/mposY/@dY",dy);
	vector<double>tof_top2;
	Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_top2']/@X_Y_Z",tof_top2);  
	for (int k=num_bars;k<num_bars+num_narrow_bars2;k++){
  	y_tof.push_back(y0+tof_top2[1]+dy*double(k-num_bars));
  	y_widths.push_back(dy);
	}
	num_bars+=num_narrow_bars2;

	// Last 19 long bars
	Get("//composition[@name='forwardTOF_top1']/mposY/@Y0",y0);
	Get("//composition[@name='forwardTOF_top1']/mposY/@dY",dy);
	vector<double>tof_top1;
	Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_top1']/@X_Y_Z",tof_top1);  
	for (int k=num_bars;k<num_bars+num_bars2;k++){
 	 	y_tof.push_back(y0+tof_top1[1]+dy*double(k-num_bars));
  	y_widths.push_back(dy);
	}
	num_bars+=num_bars2;

/*
	// two more short wide bars - IGNORE FOR NOW, ASSUME SAME Y AS OTHER SINGLE ENDED
	Get("//composition[@name='forwardTOF_south']/mposY/@Y0",y0);
	Get("//composition[@name='forwardTOF_south']/mposY/@dY",dy);
	vector<double>tof_south;
	Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_south']/@X_Y_Z",tof_south);  
	for (unsigned int k=45;k<47;k++){
  	y_tof.push_back(y0+tof_south[1]+dy*double(k-45));
	}
  */
  
return true;
2180}

2182//---------------------------------
2183// GetTargetZ
2184//---------------------------------
2185bool DGeometry::GetTargetZ(double &z_target) const
2186{
 // Default to nominal center of GlueX target
 z_target=65.;

 jgeom->SetVerbose(0);   // don't print error messages for optional detector elements

 // Check GlueX target is defined
 bool gluex_target_exists = true;
 vector<double> xyz_vessel;
 vector<double> xyz_target;
 vector<double> xyz_detector;
 if(gluex_target_exists) gluex_target_exists = Get("//composition[@name='targetVessel']/posXYZ[@volume='targetTube']/@X_Y_Z", xyz_vessel);
 if(gluex_target_exists) gluex_target_exists = Get("//composition[@name='Target']/posXYZ[@volume='targetVessel']/@X_Y_Z", xyz_target);
 if(gluex_target_exists) gluex_target_exists = Get("//posXYZ[@volume='Target']/@X_Y_Z", xyz_detector);
 if(gluex_target_exists) {
    z_target = xyz_vessel[2] + xyz_target[2] + xyz_detector[2];
 	  jgeom->SetVerbose(1);   // reenable error messages
 return true;
 }

 // Check if CPP target is defined
 bool cpp_target_exists = true;
 vector<double> xyz_TGT0;
 vector<double> xyz_TARG;
 vector<double> xyz_TargetCPP;
 if(cpp_target_exists) cpp_target_exists = Get("//composition/posXYZ[@volume='TGT0']/@X_Y_Z", xyz_TGT0);
 if(cpp_target_exists) cpp_target_exists = Get("//composition/posXYZ[@volume='TARG']/@X_Y_Z", xyz_TARG);
 if(cpp_target_exists) cpp_target_exists = Get("//composition/posXYZ[@volume='TargetCPP']/@X_Y_Z", xyz_TargetCPP);
 if(cpp_target_exists) {
    z_target = xyz_TGT0[2] + xyz_TARG[2] + xyz_TargetCPP[2];
 	  jgeom->SetVerbose(1);   // reenable error messages
    return true;
 }
 
 // Check if PrimEx Be target is defined
 bool primex_target_exists = true;
 vector<double> xyz_BETG;
 if(primex_target_exists) primex_target_exists = Get("//composition[@name='targetVessel']/posXYZ[@volume='BETG']/@X_Y_Z", xyz_BETG);
 if(primex_target_exists) primex_target_exists = Get("//composition[@name='Target']/posXYZ[@volume='targetVessel']/@X_Y_Z", xyz_target);
 if(primex_target_exists) primex_target_exists = Get("//posXYZ[@volume='Target']/@X_Y_Z", xyz_detector);
 if(primex_target_exists) {

   z_target = xyz_BETG[2] + xyz_target[2] + xyz_detector[2];

   //     cout << " PrimEx Be targer selected. Z target =   = " << z_target << endl;
   
   jgeom->SetVerbose(1);   // reenable error messages
   return true;
 }
 


 jout << " WARNING: Unable to get target location from XML for any of GlueX, PrimEx, or CPP targets. It's likely an empty target run. Using default of " << 
   z_target << " cm" << endl;

 jgeom->SetVerbose(1);   // reenable error messages

 return false;
2244}

2246//---------------------------------
2247// GetTargetLength
2248//---------------------------------
2249bool DGeometry::GetTargetLength(double &target_length) const
2250{
 jgeom->SetVerbose(0);   // don't print error messages for optional detector elements
 vector<double> zLength;
 bool good = Get("//section[@name='Target']/pcon[@name='LIH2']/real[@name='length']/[@value]", zLength);
 jgeom->SetVerbose(1);   // reenable error messages

 target_length = good ? zLength[0]:0.0;

 return good;
2259}

2261// Get vectors of positions and norm vectors for start counter from XML
2262bool DGeometry::GetStartCounterGeom(vector<vector<DVector3> >&pos,
		    vector<vector<DVector3> >&norm
		    ) const
2265{
		    
JCalibration *jcalib = dapp->GetJCalibration(runnumber);

// Check if Start Counter geometry is present
vector<double> sc_origin;
bool got_sc = Get("//posXYZ[@volume='StartCntr']/@X_Y_Z", sc_origin);
if(got_sc){
  // Offset from beam center
  vector<double>sc_origin_delta;
  Get("//posXYZ[@volume='startCntr']/@X_Y_Z", sc_origin_delta);
  double dx=sc_origin_delta[0];
  double dy=sc_origin_delta[1];

  // z-position at upstream face of scintillators.
  double z0=sc_origin[2];
  
  // Get rotation angles
  vector<double>sc_rot_angles;
  Get("//posXYZ[@volume='startCntr']/@rot", sc_rot_angles);
  double ThetaX=sc_rot_angles[0]*M_PI3.14159265358979323846/180.;
  double ThetaY=sc_rot_angles[1]*M_PI3.14159265358979323846/180.;
  Get("//posXYZ[@volume='StartCntr']/@rot", sc_rot_angles);
  //double ThetaX=sc_rot_angles[0]*M_PI/180.;
  //double ThetaY=sc_rot_angles[1]*M_PI/180.;
  double ThetaZ=sc_rot_angles[2]*M_PI3.14159265358979323846/180.;

// Get overall alignment shifts from CCDB
  map<string,double> sc_global_offsets;
  if (jcalib->Get("START_COUNTER/global_alignment_parms",sc_global_offsets)==false) {
 		// translations
 		dx += sc_global_offsets["SC_ALIGN_X"];
 		dy += sc_global_offsets["SC_ALIGN_Y"];
 		z0 += sc_global_offsets["SC_ALIGN_Z"];

// rotations
ThetaX += sc_global_offsets["SC_ALIGN_ROTX"]*M_PI3.14159265358979323846/180.;
ThetaY += sc_global_offsets["SC_ALIGN_ROTY"]*M_PI3.14159265358979323846/180.;
ThetaZ += sc_global_offsets["SC_ALIGN_ROTZ"]*M_PI3.14159265358979323846/180.;
  }

  double num_paddles;
  Get("//mposPhi[@volume='STRC']/@ncopy",num_paddles); 
  double dSCdphi = M_TWO_PI6.28318530717958647692/num_paddles;
  
  vector<vector<double> > sc_rioz;
  GetMultiple("//pgon[@name='STRC']/polyplane/@Rio_Z", sc_rioz);
  
  // Get individual paddle alignment parameters
  vector< map<string,double> > sc_paddle_offsets;
bool loaded_paddle_offsets = false;
  if (jcalib->Get("START_COUNTER/paddle_alignment_parms",sc_paddle_offsets)==false)
  	loaded_paddle_offsets = true;
  
  // Create vectors of positions and normal vectors for each paddle
  for (unsigned int i=0;i<30;i++){
    double phi=ThetaZ+dSCdphi*(double(i)+0.5);
    double sinphi=sin(phi);
    double cosphi=cos(phi);
    double r=0.5*(sc_rioz[0][0]+sc_rioz[0][1]);
    DVector3 oldray;
    // Rotate by phi and take into account the tilt
    double x=r*cosphi;//+dx;
    double y=r*sinphi;//+dy;
    DVector3 ray(x,y,sc_rioz[0][2]);
    ray.RotateX(ThetaX);
    ray.RotateY(ThetaY);
    
    // Create stl-vectors to store positions and norm vectors
    vector<DVector3>posvec;
    vector<DVector3>dirvec;
    // Loop over radial/z positions describing start counter geometry from xml
    for(unsigned int k = 1; k < sc_rioz.size(); ++k){
oldray=ray;
r=0.5*(sc_rioz[k][0]+sc_rioz[k][1]);
// Point in plane of scintillator
x=r*cosphi;//+dx;
y=r*sinphi;//+dy;
ray.SetXYZ(x,y,sc_rioz[k][2]); 
ray.RotateX(ThetaX);
ray.RotateY(ThetaY);  
// Apply alignment parameters
if(loaded_paddle_offsets) {
    // allow for a maximum extent of the paddle in z
 double max_z = sc_paddle_offsets[i]["SC_MAX_Z"];
 if(ray.Z() > max_z) {
ray.SetZ(max_z);
 }
 // allow for a modification of the bend angle of the paddle (18.5 deg from horizontal)
 // this should just be a perturbation around this angle, so assume a linear interpolation
 double delta_theta = sc_paddle_offsets[i]["SC_CURVE_THETA"];   // in degrees, r of curvature = 120 cm
 ray.SetX(ray.X()+delta_theta*1.65);   // 1 degree ~ 1.65 cm in x
 ray.SetY(ray.Y()-delta_theta*0.55);   // 1 degree ~ 0.55 cm in y
}
// Second point in the plane of the scintillator
DVector3 ray2(r,10.,sc_rioz[k][2]);
ray2.RotateZ(phi+0.5*dSCdphi*(1.+1./15.*((i>14)?29-i:i)));	
ray2.RotateX(ThetaX);
ray2.RotateY(ThetaY);  
// Compute normal vector to plane
DVector3 dir=(ray-oldray).Cross(ray2-oldray);
dir.SetMag(1.);
dirvec.push_back(dir);
posvec.push_back(DVector3(oldray.X()+dx,oldray.Y()+dy,oldray.Z()+z0));
    }
    posvec.push_back(DVector3(ray.X(),ray.Y(),ray.Z()+z0)); //SAVE THE ENDPOINT OF THE LAST PLANE
    pos.push_back(posvec);
    norm.push_back(dirvec);
  
    posvec.clear();
    dirvec.clear();
  }
  
}
return got_sc;
2380}


←

libraries/HDGEOMETRY/DGeometry.h

→

1// $Id$
2//
3//    File: DGeometry.h
4// Created: Thu Apr  3 08:43:06 EDT 2008
5// Creator: davidl (on Darwin swire-d95.jlab.org 8.11.1 i386)
6//

8#ifndef _DGeometry_
9#define _DGeometry_

11#include <pthread.h>
12#include <map>

14#include <JANA/jerror.h>
15#include <JANA/JGeometry.h>
16using namespace jana;

18#include <DANA/DApplication.h>
19#include "FDC/DFDCGeometry.h"
20#include "FDC/DFDCWire.h"
21#include "FDC/DFDCCathode.h"
22#include "CDC/DCDCWire.h"

24#include <particleType.h>
25#include <DVector3.h>
26#include "DMaterial.h"
27#include "DMaterialMap.h"
28using namespace jana;

30class DApplication;
31class DMagneticFieldMap;
32class DLorentzDeflections;


35class DGeometry{
public:
DGeometry(JGeometry *jgeom, DApplication *dapp, int32_t runnumber);
virtual ~DGeometry();
virtual const char* className(void){return static_className();}
static const char* static_className(void){return "DGeometry";}

JGeometry* GetJGeometry(void) {return jgeom;}
DMagneticFieldMap* GetBfield(void) const;
DLorentzDeflections *GetLorentzDeflections(void);

// These methods just map to the same methods in JGeometry. Ideally, we'd
// base DGeometry on JGeometry and so we'd get these automatically.
// However, that would require a more complicated generator mechanism
// where the geometry objects are made outside of JANA.
bool Get(string xpath, string &sval) const {return jgeom->Get(xpath, sval);}
bool Get(string xpath, map<string, string> &svals) const {return jgeom->Get(xpath, svals);}
template<class T> bool Get(string xpath, T &val) const {return jgeom->Get(xpath, val);}
3
←
Calling 'JGeometry::Get'→
7
←
Returning from 'JGeometry::Get'→
8
←
Returning without writing to 'val'→
template<class T> bool Get(string xpath, vector<T> &vals, string delimiter=" ") const {return jgeom->Get(xpath, vals, delimiter);}
template<class T> bool Get(string xpath, map<string,T> &vals) const {return jgeom->Get(xpath, vals);}

// The GNU 3.2.3 compiler has a problem resolving the ambiguity between
// Get(string, T&val) and Get(string, vector<T> &vals, string) above.
// This does not seem to be a problem with the 4.0 compiler. To get
// around this, some non-templated versions are provided (eeech!).
bool Get(string xpath, vector<double> &vals, string delimiter=" ") const {return jgeom->Get(xpath, vals, delimiter);}
bool Get(string xpath, vector<int> &vals, string delimiter=" ") const {return jgeom->Get(xpath, vals, delimiter);}
bool Get(string xpath, vector<float> &vals, string delimiter=" ") const {return jgeom->Get(xpath, vals, delimiter);}

bool GetMultiple(string xpath,vector<vector<double> >&vals,
		 string delimiter=" ") const
{return jgeom->GetMultiple(xpath,vals,delimiter);}

typedef struct{
  double du,dphi,dz;
    }fdc_wire_offset_t;
    typedef struct{
       double dPhiX,dPhiY,dPhiZ;
    }fdc_wire_rotation_t;
    typedef struct{
       double du,dphi;
    }fdc_cathode_offset_t;
    typedef struct{
       double dx_u,dy_u,dx_d,dy_d;
    }cdc_offset_t;

    typedef pair<string, map<string,string> > node_t;
    typedef vector<node_t> xpathparsed_t;

    void FindNodes(string xpath, vector<xpathparsed_t> &matched_xpaths) const;

    // Methods for accessing material map tables obtained from calibDB
    jerror_t FindMat(DVector3 &pos, double &rhoZ_overA, double &rhoZ_overA_logI, double &RadLen) const;
    jerror_t FindMat(DVector3 &pos, double &density, double &A, double &Z, double &RadLen) const;
    jerror_t FindMatALT1(DVector3 &pos, DVector3 &mom, double &KrhoZ_overA, 
          double &rhoZ_overA,double &LnI,
          double &X0, double *s_to_boundary=NULL__null) const;
    jerror_t FindMatKalman(const DVector3 &pos,const DVector3 &mom,
          double &KrhoZ_overA,
          double &rhoZ_overA,double &LnI,double &Z,
          double &chi2c_factor,
          double &chi2a_factor,
          double &chi2a_factor2,
          unsigned int &last_index,
          double *s_to_boundary=NULL__null) const; 
    jerror_t FindMatKalman(const DVector3 &pos,const DVector3 &mom,
	     double &KrhoZ_overA,
	     double &rhoZ_overA,double &LnI,double &Z,
	     unsigned int &last_index,
	     double *s_to_boundary=NULL__null) const;
    jerror_t FindMatKalman(const DVector3 &pos,
          double &KrhoZ_overA,
          double &rhoZ_overA,double &LnI,
          double &Z,
          double &chi2c_factor,
          double &chi2a_factor,
          double &chi2a_factor2,
          unsigned int &last_index) const; 
    jerror_t FindMatKalman(const DVector3 &pos,
          double &KrhoZ_overA,
          double &rhoZ_overA,double &LnI,
          double &Z,unsigned int &last_index) const;

    const DMaterialMap::MaterialNode* FindMatNode(DVector3 &pos) const;
    const DMaterialMap* FindDMaterialMap(DVector3 &pos) const;

    // Convenience methods
    const DMaterial* GetDMaterial(string name) const;

    bool GetFDCWires(vector<vector<DFDCWire *> >&fdcwires) const;
    bool GetFDCCathodes(vector<vector<DFDCCathode *> >&fdccathodes) const;
    bool GetFDCZ(vector<double> &z_wires) const; ///< z-locations for each of the FDC wire planes in cm
    bool GetFDCStereo(vector<double> &stereo_angles) const; ///< stereo angles of each of the FDC wire layers
    bool GetFDCRmin(vector<double> &rmin_packages) const; ///< beam hole size for each FDC package in cm
    bool GetFDCRmax(double &rmax_active_fdc) const; ///< outer radius of FDC active area in cm

    bool GetCDCWires(vector<vector<DCDCWire *> >&cdcwires) const;
    bool GetCDCOption(string &cdc_option) const; ///< get the centralDC_option-X string
    bool GetCDCCenterZ(double &cdc_center_z) const; ///< z-location of center of CDC wires in cm
    bool GetCDCAxialLength(double &cdc_axial_length) const; ///< length of CDC axial wires in cm
    bool GetCDCStereo(vector<double> &cdc_stereo) const; ///< stereo angle for each CDC layer in degrees
    bool GetCDCRmid(vector<double> &cdc_rmid) const; ///< Distance of the center of CDC wire from beamline for each layer in cm
    bool GetCDCNwires(vector<int> &cdc_nwires) const; ///< Number of wires for each CDC layer
    bool GetCDCEndplate(double &z,double &dz,double &rmin,double &rmax) const; 
    bool GetCDCAxialWires(unsigned int ring,unsigned int ncopy,
          double zcenter,double dz,
          vector<vector<cdc_offset_t> >&cdc_offsets,
          vector<DCDCWire*> &axialwires,
          vector<double>&rot_angles,double dx,
          double dy) const;
    bool GetCDCStereoWires(unsigned int ring,unsigned int ncopy,
          double zcenter,
          double dz,
          vector<vector<cdc_offset_t> >&cdc_offsets,
          vector<DCDCWire*> &stereowires,
          vector<double>&rot_angles,
          double dx,double dy) const;

    bool GetBCALRmin(float &bcal_rmin) const; ///< minimum distance of BCAL module from beam line
    bool GetBCALfADCRadii(vector<float> &fADC_radii) const; ///< fADC radii including the outer radius of the last layer
    bool GetBCALNmodules(unsigned int &bcal_nmodules) const; ///< Number of BCAL modules
    bool GetBCALCenterZ(float &bcal_center_z) const; ///< z-location of center of BCAL module in cm
    bool GetBCALLength(float &bcal_length) const; ///< length of BCAL module in cm
    bool GetBCALDepth(float &bcal_depth) const; ///< depth (or height) of BCAL module in cm
    bool GetBCALPhiShift(float &bcal_phi_shift) const; ///< phi angle in degrees that first BCAL module is shifted from being centered at ph=0.0

    bool GetCCALZ(double &z_ccal) const; /// z-location of front face of CCAL in cm

    bool GetFCALZ(double &z_fcal) const; ///< z-location of front face of FCAL in cm
    bool GetDIRCZ(double &z_dirc) const; ///< z-location of DIRC in cm
    bool GetTOFZ(vector<double> &z_tof) const; ///< z-location of front face of each of TOF in cm
 bool GetTOFPaddlePerpPositions(vector<double> &y_tof, vector<double> &y_widths) const;
 bool GetTOFPaddleParameters(map<string,double> &paddle_params) const;
    bool GetTargetZ(double &z_target) const; ///< z-location of center of target
    bool GetTargetLength(double &target_length) const; ///< z-location of center of target

    bool GetTRDZ(vector<double> &z_trd) const; ///< z-locations for each of the TRD/GEM planes in cm
    
    bool GetFCALPosition(double &x,double &y,double &z) const;
    bool GetCCALPosition(double &x,double &y,double &z) const;

    bool GetFCALInsertRowSize(int &insert_row_size) const;
    bool GetFCALBlockSize(vector<double> &block) const;
    bool GetFCALInsertBlockSize(vector<double> &block) const;

    bool GetStartCounterGeom(vector<vector<DVector3> >&pos,
          vector<vector<DVector3> >&norm) const; // < vectors containing positions and norm 3-vectors for start counter 
    // There are 30 sets of positions (pos) of points along the 
    // start counter, one set for each paddle, and a corresponding
    // set of norms at various points along the start counter.
    // For example, to access the most upstream point of paddle 1
    // use pos[0][0].  The end of the barrel section before the 
    // nose region is at pos[0][1].  The tip of the nose region 
    // for this paddle is at pos[0][pos[0].size()-1].  The bend
    // region is modeled by many closely-spaced points starting 
    // after pos[0][1].

    vector<DMaterialMap*> GetMaterialMapVector(void) const;

 protected:
    DGeometry(){}
    void ReadMaterialMaps(void) const;
    void GetMaterials(void) const;
    bool GetCompositeMaterial(const string &name, double &density, double &radlen) const;

 private:
    JGeometry *jgeom;
    DApplication *dapp;
    mutable DMagneticFieldMap *bfield;
    int32_t runnumber;
    mutable vector<DMaterial*> materials;			/// Older implementation to keep track of material specs without ranges
    mutable vector<DMaterialMap*> materialmaps;	/// Material maps generated automatically(indirectly) from XML with ranges and specs
    mutable bool materialmaps_read;
    mutable bool materials_read;

    mutable pthread_mutex_t bfield_mutex;
    mutable pthread_mutex_t materialmap_mutex;
    mutable pthread_mutex_t materials_mutex;


215};

217#endif // _DGeometry_


←

/w/halld-scifs17exp/halld2/home/sdobbs/Software/jana/jana_0.8.2/Linux_CentOS7.7-x86_64-gcc4.8.5/include/JANA/JGeometry.h

1// $Id$
2//
3//    File: JGeometry.h
4// Created: Fri Jul  6 16:24:24 EDT 2007
5// Creator: davidl (on Darwin fwing-dhcp61.jlab.org 8.10.1 i386)
6//

8#ifndef _JGeometry_
9#define _JGeometry_

11#include "jerror.h"

13#include <map>
14#include <string>
15#include <sstream>
16#include <vector>
17using std::map;
18using std::string;
19using std::stringstream;
20using std::vector;

22// Place everything in JANA namespace
23namespace jana{

25/// JGeometry is a virtual base class used to define the interface by
26/// which geometry information can be obtained in JANA.
27/// Implementing this base class allows the JANA end user to be
28/// agnostic as to the details of how the geometry info is stored.
29/// The geometry can be stored in a database or any number of file
30/// formats. The files can be stored locally, or on the network
31/// somewhere.
32///
33/// The primary advantage here is that it allows one to work with
34/// local files, but then easily switch to a remote source method
35/// when appropriate without requiring modifications to the end
36/// user code.
37///
38/// On the user side they will call one of the Get(...) methods which all get
39/// translated into a call of one of the two vitural methods:
40///
41///   virtual bool Get(string namepath, string &sval, map<string, string> &where)=0;
42///   virtual bool Get(string namepath, map<string, string> &svals, map<string, string> &where)=0;
43///
44/// These two virtual methods along with one to get a list of the available
45/// namepaths are the only things that need to be implemented
46/// in a concrete subclass of JGeometry.
47///
48/// A geometry element is specified by its <i>namepath</i> and an optional
49/// set of qualifiers (the <i>where</i> argument). The <i>namepath</i>
50/// is a hierarchal list of elements separated by forward slashes(/)
51/// analogous to a path on a unix filesystem. This path is always assumed
52/// to be relative to the <i>url</i> specified in the constructor. So,
53/// for instance, suppose one kept the geometry in a set XML files on the
54/// local filesystem and wished to access information from the file
55///
56/// <tt>/home/joe/calib/geom_Oct10_2017.xml</tt>
57///
58/// One would specify the <i>url</i> as:
59///
60/// file:///home/joe/calib/geom_Oct10_2017.xml
61///
62/// and then the namepath could be specified as the string:
63///
64/// "TOF/bar/X_Y_Z"
65///
66/// which would indicate the attribute <i>"X_Y_Z"</i> of the
67/// subtag <i>"bar"</i> of the tag <i>"TOF"</i> in the file
68/// <i>"/home/joe/calib/geom_Oct10_2017.xml"</i>

70class JGeometry{
public:
JGeometry(string url, int run, string context="default"){
	this->url = url;
	this->run_requested = run;
	this->context = context;
}
virtual ~JGeometry(){}
virtual const char* className(void){return static_className();}
static const char* static_className(void){return "JGeometry";}

typedef enum{
	// Used to specify which (if any) attributes should be included in
	// the values obtained through GetXPaths().
	attr_level_none = 0,	// Don't include any attributes. Only node names.
	attr_level_last = 1,	// Include the attributes for the last node only.
	attr_level_all  = 2	// Include attributes for all nodes.
}ATTR_LEVEL_t;

void SetVerbose(int newval){ verbose = newval; }

// Virtual methods called through base class
virtual bool Get(string xpath, string &sval)=0;
virtual bool Get(string xpath, map<string, string> &svals)=0;
virtual bool GetMultiple(string xpath, vector<string> &vsval)=0;
virtual bool GetMultiple(string xpath, vector<map<string, string> >&vsvals)=0;
virtual void GetXPaths(vector<string> &xpaths, ATTR_LEVEL_t level=attr_level_last, const string &filter="")=0;
virtual string GetChecksum(void) const {return string("not supported");}

// Templated methods that can return more useful forms
template<class T> bool Get(string xpath, T &val);
template<class T> bool Get(string xpath, vector<T> &vals, string delimiter=" ");
template<class T> bool Get(string xpath, map<string,T> &vals);
template<class T> bool GetMultiple(string xpath, vector<T> &vval);
template<class T> bool GetMultiple(string xpath, vector<vector<T> > &vvals, string delimiter=" ");
template<class T> bool GetMultiple(string xpath, vector<map<string,T> > &vvals);

const int& GetRunRequested(void) const {return run_requested;}
const int& GetRunFound(void) const {return run_found;}
const int& GetRunMin(void) const {return run_min;}
const int& GetRunMax(void) const {return run_max;}
const string& GetContext(void) const {return context;}
const string& GetURL(void) const {return url;}

protected:
int run_min;
int run_max;
int run_found;

int verbose=1;

private:
JGeometry(){} // Don't allow trivial constructor

int run_requested;
string context;
string url;
map<string,string> anywhere; // an empty map means no constraints

129};


132//-------------
133// Get  (single version)
134//-------------
135template<class T>
136bool JGeometry::Get(string xpath, T &val)
137{
/// Templated method used to get a single geometry element.
///
/// This method will get the specified geometry element in the form of
/// a string using the virtual (non-templated) Get(...) method. It will
/// then convert the string into the data type on which <i>val</i> is
/// based. It does this using the stringstream
/// class so T is restricted to the types stringstream understands (int, float, 
/// double, string, ...).
///
/// If no element of the specified name is found, a value
/// of boolean "false" is returned. A value of "true" is 
/// returned upon success.

// Get values in the form of a string
string sval;
bool res = Get(xpath, sval);
if(!res)return res;
4
←
Assuming 'res' is false→
5
←
Taking true branch→
6
←
Returning without writing to 'val'→

// Convert the string to type "T" and copy it into val.
// Use stringstream to convert from a string to type "T"
stringstream ss(sval);
ss >> val;

return res;
162}

164//-------------
165// Get  (vector version)
166//-------------
167template<class T>
168bool JGeometry::Get(string xpath, vector<T> &vals, string delimiter)
169{
/// Templated method used to get a set of values from a geometry attribute.
///
/// This method can be used to get a list of values (possibly only one)
/// from a single geometry attribute specified by xpath. The attribute
/// is obtained as a string using the non-templated Get(...) method
/// and the string broken into tokens separated by the delimiter
/// (which defaults to a single white space). Each
/// token is then converted into type T using the stringstream class
/// so T is restricted to the types stringstream understands (int, float, 
/// double, string, ...).
///
/// If no element of the specified name is found, a value
/// of boolean "false" is returned. A value of "true" is 
/// returned upon success.

// Get values in the form of strings
vals.clear();
string svals;
bool res = Get(xpath, svals);
if(!res)return res;

string::size_type pos_start = svals.find_first_not_of(delimiter,0);
while(pos_start != string::npos){
string::size_type pos_end = svals.find_first_of(delimiter, pos_start);
if(pos_end==string::npos)pos_end=svals.size();

T v;
string val = svals.substr(pos_start, pos_end-pos_start);
stringstream ss(val);
ss >> v;
vals.push_back(v);

pos_start = svals.find_first_not_of(delimiter, pos_end);
}

return res;
206}

208//-------------
209// Get  (map version)
210//-------------
211template<class T>
212bool JGeometry::Get(string xpath, map<string,T> &vals)
213{
/// Templated method used to get a set of geometry attributes.
///
/// This method can be used to get a list of all attributes for
/// a given xpath. The attributes are copied into the <i>vals</i>
/// map with the attribute name as the key and the attribute
/// value as the value. This relies on the non-templated, virtual
/// Get(string, map<string,string>&) method to first get the values
/// in the form of strings. It converts them using the stringstream
/// class so T is restricted to the types it understands (int, float, 
/// double, string, ...).
///
/// If no element of the specified name is found, a value
/// of boolean "false" is returned. A value of "true" is 
/// returned upon success.

// Get values in the form of strings
map<string, string> svals;
bool res = Get(xpath, svals);

// Loop over values, converting the strings to type "T" and
// copying them into the vals map.
vals.clear();
map<string,string>::const_iterator iter;
for(iter=svals.begin(); iter!=svals.end(); ++iter){
// Use stringstream to convert from a string to type "T"
T v;
stringstream ss(iter->second);
ss >> v;
vals[iter->first] = v;
}

return res;
246}


249//-------------
250// GetMultiple  (single version)
251//-------------
252template<class T>
253bool JGeometry::GetMultiple(string xpath, vector<T> &vval)
254{
/// Templated method used to get multiple entries satisfying a single xpath.
///
/// This method will get the specified geometry element in the form of
/// a string using the virtual (non-templated) Get(...) method. It will
/// then convert the string into the data type on which <i>val</i> is
/// based. It does this using the stringstream
/// class so T is restricted to the types stringstream understands (int, float, 
/// double, string, ...).
///
/// This differs from the similar Get() method in that the geometry tree
/// will be searched for all nodes satisfying the given xpath and all
/// all values will be copied into the container provided. In Get(), only
/// the first node encountered that satisfies the xpath will be copied.
///
/// If no element of the specified name is found, a value
/// of boolean "false" is returned. A value of "true" is 
/// returned upon success.

// Get values in the form of a string
vector<string> vsval;
bool res = GetMultiple(xpath, vsval);
if(!res)return res;

// Convert the string to type "T" and copy it into val.
// Use stringstream to convert from a string to type "T"
for(unsigned int i=0; i<vsval.size(); i++){
stringstream ss(vsval[i]);
T val;
ss >> val;
vval.push_back(val);
}

return res;
288}

290//-------------
291// GetMultiple  (vector version)
292//-------------
293template<class T>
294bool JGeometry::GetMultiple(string xpath, vector<vector<T> > &vvals, string delimiter)
295{
/// Templated method used to get a set of values from a geometry attribute.
///
/// This method can be used to get a list of values (possibly only one)
/// from a single geometry attribute specified by xpath. The attribute
/// is obtained as a string using the non-templated Get(...) method
/// and the string broken into tokens separated by the delimiter
/// (which defaults to a single white space). Each
/// token is then converted into type T using the stringstream class
/// so T is restricted to the types stringstream understands (int, float, 
/// double, string, ...).
///
/// If no element of the specified name is found, a value
/// of boolean "false" is returned. A value of "true" is 
/// returned upon success.

// Get values in the form of strings
vvals.clear();
vector<string> vsvals;
bool res = GetMultiple(xpath, vsvals);
if(!res)return res;

for(unsigned int i=0; i<vsvals.size(); i++){
string &svals = vsvals[i];

string::size_type pos_start = svals.find_first_not_of(delimiter,0);
vector<T> vals;
while(pos_start != string::npos){
	string::size_type pos_end = svals.find_first_of(delimiter, pos_start);
	if(pos_end==string::npos)pos_end=svals.size();

	T v;
	string val = svals.substr(pos_start, pos_end-pos_start);
	stringstream ss(val);
	ss >> v;
	vals.push_back(v);
	
	pos_start = svals.find_first_not_of(delimiter, pos_end);
}

vvals.push_back(vals);
}

return res;
339}

341//-------------
342// GetMultiple  (map version)
343//-------------
344template<class T>
345bool JGeometry::GetMultiple(string xpath, vector<map<string,T> > &vvals)
346{
/// Templated method used to get a set of geometry attributes.
///
/// This method can be used to get a list of all attributes for
/// a given xpath. The attributes are copied into the <i>vals</i>
/// map with the attribute name as the key and the attribute
/// value as the value. This relies on the non-templated, virtual
/// Get(string, map<string,string>&) method to first get the values
/// in the form of strings. It converts them using the stringstream
/// class so T is restricted to the types it understands (int, float, 
/// double, string, ...).
///
/// If no element of the specified name is found, a value
/// of boolean "false" is returned. A value of "true" is 
/// returned upon success.

// Get values in the form of strings
vector<map<string, string> > vsvals;
bool res = GetMultiple(xpath, vsvals);

// Loop over values, converting the strings to type "T" and
// copying them into the vals map.
vvals.clear();

for(unsigned int i=0; i<vsvals.size(); i++){
map<string,string> &svals = vsvals[i];

map<string,T> vals;
map<string,string>::const_iterator iter;
for(iter=svals.begin(); iter!=svals.end(); ++iter){
	// Use stringstream to convert from a string to type "T"
	T v;
	stringstream ss(iter->second);
	ss >> v;
	vals[iter->first] = v;
}

vvals.push_back(vals);
}

return res;
387}

389} // Close JANA namespace

391#endif // _JGeometry_