Hall-D Software  alpha
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
DGeometry.cc
Go to the documentation of this file.
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 //
7 
8 #include <algorithm>
9 using namespace std;
10 
11 #include "DGeometry.h"
12 #include <JANA/JGeometryXML.h>
13 #include "FDC/DFDCWire.h"
14 #include "FDC/DFDCGeometry.h"
15 #include <ansi_escape.h>
16 
17 using namespace std;
18 
19 #ifndef M_TWO_PI
20 #define M_TWO_PI 6.28318530717958647692
21 #endif
22 
23 //---------------------------------
24 // DGeometry (Constructor)
25 //---------------------------------
26 DGeometry::DGeometry(JGeometry *jgeom, DApplication *dapp, int32_t runnumber)
27 {
28  this->jgeom = jgeom;
29  this->dapp = dapp;
30  this->bfield = NULL; // don't ask for B-field object before we're asked for it
31  this->runnumber = runnumber;
32  this->materialmaps_read = false;
33  this->materials_read = false;
34 
35  pthread_mutex_init(&bfield_mutex, NULL);
36  pthread_mutex_init(&materialmap_mutex, NULL);
37  pthread_mutex_init(&materials_mutex, NULL);
38 
39  ReadMaterialMaps();
40 }
41 
42 //---------------------------------
43 // ~DGeometry (Destructor)
44 //---------------------------------
46 {
47  pthread_mutex_lock(&materials_mutex);
48  for(unsigned int i=0; i<materials.size(); i++)delete materials[i];
49  materials.clear();
50  pthread_mutex_unlock(&materials_mutex);
51 
52  pthread_mutex_lock(&materialmap_mutex);
53  for(unsigned int i=0; i<materialmaps.size(); i++)delete materialmaps[i];
54  materialmaps.clear();
55  pthread_mutex_unlock(&materialmap_mutex);
56 }
57 
58 //---------------------------------
59 // GetBfield
60 //---------------------------------
62 {
63  pthread_mutex_lock(&bfield_mutex);
64  if(bfield == NULL) bfield = dapp->GetBfield(runnumber);
65  pthread_mutex_unlock(&bfield_mutex);
66 
67  return bfield;
68 }
69 
70 //---------------------------------
71 // GetLorentzDeflections
72 //---------------------------------
74 {
75  return dapp->GetLorentzDeflections(runnumber);
76 }
77 
78 //---------------------------------
79 // GetMaterialMapVector
80 //---------------------------------
81 vector<DMaterialMap*> DGeometry::GetMaterialMapVector(void) const
82 {
83 // ReadMaterialMaps();
84 
85  return materialmaps;
86 }
87 
88 //---------------------------------
89 // ReadMaterialMaps
90 //---------------------------------
92 {
93  /// This gets called by several "FindMat" methods below. It
94  /// will return immediately if the material maps have already
95  /// been read in. If they have not been read in, then it reads
96  /// them and sets a flag so that they are only read in once.
97  ///
98  /// Orginally, this code resided in the constructor, but was
99  /// moved here so that programs not requiring the material
100  /// maps could start up quicker and skip reading them in altogether.
101 
102  // Lock mutex so we can check flag to see if maps have already
103  // been read in
104  pthread_mutex_lock(&materialmap_mutex);
105  if(materialmaps_read){
106  // Maps are already read. Unlock mutex and return now
107  pthread_mutex_unlock(&materialmap_mutex);
108  return;
109  }
110 
111  JCalibration * jcalib = dapp->GetJCalibration(runnumber);
112  if(!jcalib){
113  _DBG_<<"ERROR: Unable to get JCalibration object!"<<endl;
114  pthread_mutex_unlock(&materialmap_mutex);
115  return;
116  }
117 
118  // Get a list of all namepaths so we can find all of the material maps.
119  // We want to read in all material maps with a standard naming scheme
120  // so the number and coverage of the piece-wise maps can be changed
121  // without requiring recompilation.
122  //vector<DMaterialMap*> material_maps;
123  vector<string> namepaths;
124  jcalib->GetListOfNamepaths(namepaths);
125  vector<string> material_namepaths;
126  for(unsigned int i=0; i<namepaths.size(); i++){
127  if(namepaths[i].find("Material/material_map")==0)material_namepaths.push_back(namepaths[i]);
128  }
129 
130  // Sort alphabetically so user controls search sequence by map name
131  sort(material_namepaths.begin(), material_namepaths.end());
132 
133  // Inform user what's happening
134  if(material_namepaths.size()==0){
135  jerr<<"No material maps found in calibration DB!!"<<endl;
136  pthread_mutex_unlock(&materialmap_mutex);
137  return;
138  }
139  jout<<"Found "<<material_namepaths.size()<<" material maps in calib. DB"<<endl;
140 
141  if(false){ // save this to work off configuration parameter
142  jout<<"Will read in the following:"<<endl;
143  for(unsigned int i=0; i<material_namepaths.size(); i++){
144  jout<<" "<<material_namepaths[i]<<endl;
145  }
146  }
147 
148  // Actually read in the maps
149  uint32_t Npoints_total=0;
150  //cout<<endl; // make empty line so material map can overwrite it below
151  for(unsigned int i=0; i<material_namepaths.size(); i++){
152  // DMaterialMap constructor prints line so we conserve real
153  // estate by having each recycle the line
154  //cout<<ansi_up(1)<<string(85, ' ')<<"\r";
155  DMaterialMap *mat = new DMaterialMap(material_namepaths[i], jcalib);
156  if( ! mat->IS_VALID ) {
157  // This particular map may not exist for this run/variation
158  // (e.g. CPP uses maps downstream of TOF)
159  delete mat;
160  continue;
161  }
162  materialmaps.push_back(mat);
163  Npoints_total += (unsigned int)(mat->GetNr()*mat->GetNz());
164  }
165  //cout<<ansi_up(1)<<string(85, ' ')<<"\r";
166  jout<<"Read in "<<materialmaps.size()<<" material maps for run "<<runnumber<<" containing "<<Npoints_total<<" grid points total"<<endl;
167 
168  // Set flag that maps have been read and unlock mutex
169  materialmaps_read = true;
170  pthread_mutex_unlock(&materialmap_mutex);
171 }
172 
173 //---------------------------------
174 // FindNodes
175 //---------------------------------
176 void DGeometry::FindNodes(string xpath, vector<xpathparsed_t> &matched_xpaths) const
177 {
178  /// Find all nodes that match the specified xpath and return them as
179  /// fully parsed lists of the nodes and attributes.
180  ///
181  /// The matched_xpaths variable has 4 levels of STL containers nested
182  /// together! The node_t data type contains 2 of them and represents
183  /// a single tag with the "first" member being the tag name
184  /// and the "second" member being a map of all of the attributes
185  /// of that tag along with their values.
186  ///
187  /// The xpathparsed_t data type is a STL vector of node_t objects
188  /// that comprises a complete xpath. The data type of matched_xpaths
189  /// is a STL vector if xpathparsed_t objects and so comprises a
190  /// list of complete xpaths.
191 
192  /// We do this by calling the GetXPaths() method of JGeometry to get
193  /// a list of all xpaths. Then we pass all of those in to
194  /// JGeometryXML::ParseXPath() to get a parsed list for each. This
195  /// is compared to the parsed values of the xpath passed to us
196  /// (also parsed by JGeometryXML::ParseXPath()) to find matches.
197 
198  // Make sure matched_xpaths is empty
199  matched_xpaths.clear();
200 
201  // Cast JGeometry into a JGeometryXML
202  JGeometryXML *jgeomxml = dynamic_cast<JGeometryXML*>(jgeom);
203 
204  // Parse our target xpath
205  xpathparsed_t target;
206  string unused_string;
207  unsigned int unused_int;
208  jgeomxml->ParseXPath(xpath, target, unused_string, unused_int);
209 
210  // Get all xpaths for current geometry source
211  vector<string> allxpaths;
212  jgeom->GetXPaths(allxpaths, JGeometry::attr_level_all);
213 
214  // Loop over xpaths
215  for(unsigned int i=0; i<allxpaths.size(); i++);
216 }
217 
218 //---------------------------------
219 // FindMatALT1 - parameters for alt1 fitter.
220 //---------------------------------
221 jerror_t DGeometry::FindMatALT1(DVector3 &pos, DVector3 &mom,double &KrhoZ_overA,
222  double &rhoZ_overA,double &LnI,
223  double &X0, double *s_to_boundary) const
224 {
225 // ReadMaterialMaps();
226 
227  for(unsigned int i=0; i<materialmaps.size(); i++){
228  jerror_t err = materialmaps[i]->FindMatALT1(pos,KrhoZ_overA, rhoZ_overA,LnI,X0);
229  if(err==NOERROR){
230  // We found the material map containing this point. If a non-NULL
231  // pointer was passed in for s_to_boundary, then search through all
232  // material maps above and including this one to find the estimated
233  // distance to the nearest boundary the swimmer should step to. Maps
234  // after this one would only be considered outside of this one so
235  // there is no need to check them.
236  if(s_to_boundary==NULL)return NOERROR; // User doesn't want distance to boundary
237  *s_to_boundary = 1.0E6;
238  for(unsigned int j=0; j<=i; j++){
239  double s = materialmaps[j]->EstimatedDistanceToBoundary(pos, mom);
240  if(s<*s_to_boundary)*s_to_boundary = s;
241  }
242  return NOERROR;
243  }
244  }
245  return RESOURCE_UNAVAILABLE;
246 }
247 
248 
249 
250 
251 //---------------------------------
252 // FindMatKalman - Kalman filter needs slightly different set of parms.
253 //---------------------------------
254 jerror_t DGeometry::FindMatKalman(const DVector3 &pos,const DVector3 &mom,
255  double &KrhoZ_overA,
256  double &rhoZ_overA,
257  double &LnI,double &Z,
258  double &chi2c_factor,double &chi2a_factor,
259  double &chi2a_corr,
260  unsigned int &last_index,
261  double *s_to_boundary) const
262 {
263 // ReadMaterialMaps();
264 
265  //last_index=0;
266  for(unsigned int i=last_index; i<materialmaps.size(); i++){
267  jerror_t err = materialmaps[i]->FindMatKalman(pos,KrhoZ_overA,
268  rhoZ_overA,LnI,chi2c_factor,
269  chi2a_factor,chi2a_corr,Z);
270  if(err==NOERROR){
271  if(i==materialmaps.size()-1) last_index=0;
272  else last_index=i;
273  if(s_to_boundary==NULL)return NOERROR; // User doesn't want distance to boundary
274 
275  *s_to_boundary = 1.0E6;
276  // If we are in the main mother volume, search through all the maps for
277  // the nearest boundary
278  if(last_index==0){
279  for(unsigned int j=0; j<materialmaps.size();j++){
280  double s = materialmaps[j]->EstimatedDistanceToBoundary(pos, mom);
281  if(s<*s_to_boundary){
282  *s_to_boundary = s;
283  }
284  }
285  }
286  else{
287  // otherwise, we found the material map containing this point.
288  double s = materialmaps[last_index]->EstimatedDistanceToBoundary(pos, mom);
289  if(s<*s_to_boundary)*s_to_boundary = s;
290  }
291  return NOERROR;
292  }
293  }
294 
295 return RESOURCE_UNAVAILABLE;
296 }
297 
298 //---------------------------------
299 jerror_t DGeometry::FindMatKalman(const DVector3 &pos,
300  double &KrhoZ_overA,
301  double &rhoZ_overA,
302  double &LnI, double &Z,double &chi2c_factor,
303  double &chi2a_factor, double &chi2a_corr,
304  unsigned int &last_index) const
305 {
306 // ReadMaterialMaps();
307 
308  //last_index=0;
309  for(unsigned int i=last_index; i<materialmaps.size(); i++){
310  jerror_t err = materialmaps[i]->FindMatKalman(pos,KrhoZ_overA,
311  rhoZ_overA,LnI,
312  chi2c_factor,chi2a_factor,
313  chi2a_corr,Z);
314  if(err==NOERROR){
315  if(i==materialmaps.size()-1) last_index=0;
316  else last_index=i;
317  return err;
318  }
319  }
320 
321  return RESOURCE_UNAVAILABLE;
322 }
323 
324 //---------------------------------
325 // FindMat
326 //---------------------------------
327 jerror_t DGeometry::FindMat(DVector3 &pos, double &rhoZ_overA, double &rhoZ_overA_logI, double &RadLen) const
328 {
329 // ReadMaterialMaps();
330 
331  for(unsigned int i=0; i<materialmaps.size(); i++){
332  jerror_t err = materialmaps[i]->FindMat(pos, rhoZ_overA, rhoZ_overA_logI, RadLen);
333  if(err==NOERROR)return NOERROR;
334  }
335  return RESOURCE_UNAVAILABLE;
336 }
337 
338 //---------------------------------
339 // FindMat
340 //---------------------------------
341 jerror_t DGeometry::FindMat(DVector3 &pos, double &density, double &A, double &Z, double &RadLen) const
342 {
343 // ReadMaterialMaps();
344 
345  for(unsigned int i=0; i<materialmaps.size(); i++){
346  jerror_t err = materialmaps[i]->FindMat(pos, density, A, Z, RadLen);
347  if(err==NOERROR)return NOERROR;
348  }
349  return RESOURCE_UNAVAILABLE;
350 }
351 
352 //---------------------------------
353 // FindMatNode
354 //---------------------------------
355 //const DMaterialMap::MaterialNode* DGeometry::FindMatNode(DVector3 &pos) const
356 //{
357 //
358 //}
359 
360 //---------------------------------
361 // FindDMaterialMap
362 //---------------------------------
364 {
365 // ReadMaterialMaps();
366 
367  for(unsigned int i=0; i<materialmaps.size(); i++){
368  const DMaterialMap* map = materialmaps[i];
369  if(map->IsInMap(pos))return map;
370  }
371  return NULL;
372 }
373 
374 //====================================================================
375 // Convenience Methods
376 //
377 // Below are defined some methods to make it easy to extract certain
378 // key values about the GlueX detector geometry from the XML source.
379 // Note that one can still use the generic Get(string xpath, ...)
380 // methods. This just packages some of them up for convenience.
381 //
382 // The one real gotcha here is that these methods must be kept
383 // in sync with the XML structure by hand. If volumes are renamed
384 // or their location within the hierachy is modified, then these
385 // routines will need to be modified as well. That, or course, is
386 // also true if you are using the generic Get(...) methods.
387 //
388 // What these methods are useful for is when minor changes are made
389 // to the XML (such as the locations of the FDC packages) they
390 // are automatically reflected here.
391 //====================================================================
392 
393 //---------------------------------
394 // GetDMaterial
395 //---------------------------------
396 const DMaterial* DGeometry::GetDMaterial(string name) const
397 {
398  /// Get a pointer to the DMaterial object with the specified name.
399  // Only fill materials member when one is actually requested
400  // and then, only fill it once.
401 
402  // Lock mutex so we can check flag to see if maps have already
403  // been read in
404  pthread_mutex_lock(&materials_mutex);
405  if(!materials_read) GetMaterials();
406  pthread_mutex_unlock(&materials_mutex);
407 
408  for(unsigned int i=0; i<materials.size(); i++){
409  if(materials[i]->GetName() == name)return materials[i];
410  }
411 
412  //_DBG_<<"No material \""<<name<<"\" found ("<<materials.size()<<" materials defined)"<<endl;
413 
414  return NULL;
415 }
416 
417 //---------------------------------
418 // GetMaterials
419 //---------------------------------
420 void DGeometry::GetMaterials(void) const
421 {
422  /// Read in all of the materials from the geometry source and create
423  /// a DMaterial object for each one.
424 
425  //=========== elements ===========
426  string filter = "//materials/element/real[@name=\"radlen\"]";
427 
428  // Get list of all xpaths
429  vector<string> xpaths;
430  jgeom->GetXPaths(xpaths, JGeometry::attr_level_all, filter);
431 
432  // Look for xpaths that have "/materials[" in them
433  for(unsigned int i=0; i<xpaths.size(); i++){
434  // Get start of "element" section
435  string::size_type pos = xpaths[i].find("/element[");
436  if(pos == string::npos)continue;
437 
438  // Get name attribute
439  string::size_type start_pos = xpaths[i].find("@name=", pos);
440  start_pos = xpaths[i].find("'", start_pos);
441  string::size_type end_pos = xpaths[i].find("'", start_pos+1);
442  if(end_pos==string::npos)continue;
443  string name = xpaths[i].substr(start_pos+1, end_pos-(start_pos+1));
444 
445  // Get values
446  char xpath[256];
447 
448  double A;
449  sprintf(xpath,"//materials/element[@name='%s']/[@a]", name.c_str());
450  if(!Get(xpath, A))continue;
451 
452  double Z;
453  sprintf(xpath,"//materials/element[@name='%s']/[@z]", name.c_str());
454  if(!Get(xpath, Z))continue;
455 
456  double density;
457  sprintf(xpath,"//materials/element[@name='%s']/real[@name='density']/[@value]", name.c_str());
458  if(!Get(xpath, density))continue;
459 
460  double radlen;
461  sprintf(xpath,"//materials/element[@name='%s']/real[@name='radlen']/[@value]", name.c_str());
462  if(!Get(xpath, radlen))continue;
463 
464  DMaterial *mat = new DMaterial(name, A, Z, density, radlen);
465  materials.push_back(mat);
466 
467  cout<<mat->toString();
468  }
469 
470  //=========== composites ===========
471  filter = "//materials/composite[@name]";
472 
473  // Get list of all xpaths
474  jgeom->GetXPaths(xpaths, JGeometry::attr_level_all, filter);
475 
476  // Look for xpaths that have "/materials[" in them
477  for(unsigned int i=0; i<xpaths.size(); i++){
478  // Get start of "composite" section
479  string::size_type pos = xpaths[i].find("/composite[");
480  if(pos == string::npos)continue;
481 
482  // Get name attribute
483  string::size_type start_pos = xpaths[i].find("@name=", pos);
484  start_pos = xpaths[i].find("'", start_pos);
485  string::size_type end_pos = xpaths[i].find("'", start_pos+1);
486  if(end_pos==string::npos)continue;
487  string name = xpaths[i].substr(start_pos+1, end_pos-(start_pos+1));
488 
489  if(GetDMaterial(name))continue; // skip duplicates
490 
491  // Get values
492  char xpath[256];
493 
494  // We should calculate an effective A and Z .... but we don't
495  bool found_all=true;
496  double A=0;
497  double Z=0;
498 
499  double density;
500  sprintf(xpath,"//materials/composite[@name='%s']/real[@name='density']/[@value]", name.c_str());
501  found_all &= Get(xpath, density);
502 
503  double radlen;
504  sprintf(xpath,"//materials/composite[@name='%s']/real[@name='radlen']/[@value]", name.c_str());
505  found_all &= Get(xpath, radlen);
506 
507  // If we didn't find the info we need (radlen and density) in the
508  // composite tag itself. Try calculating them from the components
509  if(!found_all)found_all = GetCompositeMaterial(name, density, radlen);
510 
511  // If we weren't able to get the values needed to make the DMaterial object
512  // then skip this one.
513  if(!found_all)continue;
514 
515  DMaterial *mat = new DMaterial(name, A, Z, density, radlen);
516  materials.push_back(mat);
517 
518  cout<<mat->toString();
519  }
520 
521  materials_read = true;
522 }
523 
524 //---------------------------------
525 // GetCompositeMaterial
526 //---------------------------------
527 bool DGeometry::GetCompositeMaterial(const string &name, double &density, double &radlen) const
528 {
529  // Get list of all xpaths with "addmaterial" and "fractionmass"
530  char filter[512];
531  sprintf(filter,"//materials/composite[@name='%s']/addmaterial/fractionmass[@fraction]", name.c_str());
532  vector<string> xpaths;
533  jgeom->GetXPaths(xpaths, JGeometry::attr_level_all, filter);
534 
535  // Loop over components of this composite
536 _DBG_<<"Components for compsite "<<name<<endl;
537  for(unsigned int i=0; i<xpaths.size(); i++){
538  // Get component material name
539  string::size_type start_pos = xpaths[i].find("@material=", 0);
540  start_pos = xpaths[i].find("'", start_pos);
541  string::size_type end_pos = xpaths[i].find("'", start_pos+1);
542  if(end_pos==string::npos)continue;
543  string mat_name = xpaths[i].substr(start_pos+1, end_pos-(start_pos+1));
544 
545  // Get component mass fraction
546  start_pos = xpaths[i].find("fractionmass[", 0);
547  start_pos = xpaths[i].find("@fraction=", start_pos);
548  start_pos = xpaths[i].find("'", start_pos);
549  end_pos = xpaths[i].find("'", start_pos+1);
550  if(end_pos==string::npos)continue;
551  string mat_frac_str = xpaths[i].substr(start_pos+1, end_pos-(start_pos+1));
552  double fractionmass = atof(mat_frac_str.c_str());
553 
554  _DBG_<<" "<<xpaths[i]<<" fractionmass="<<fractionmass<<endl;
555  }
556 
557  return true;
558 }
559 
560 //---------------------------------
561 // GetTraversedMaterialsZ
562 //---------------------------------
563 //void DGeometry::GetTraversedMaterialsZ(double q, const DVector3 &pos, const DVector3 &mom, double z_end, vector<DMaterialStep> &materialsteps)
564 //{
565  /// Find the materials traversed by a particle swimming from a specific
566  /// position with a specific momentum through the magnetic field until
567  /// it reaches a specific z-location. Energy loss is not included in
568  /// the swimming since this method itself is assumed to be one of the
569  /// primary means of determining energy loss. As such, one should not
570  /// pass in a value of z_end that is far from pos.
571  ///
572  /// The vector materialsteps will be filled with DMaterialStep objects
573  /// corresponding to each of the materials traversed.
574  ///
575  /// The real work here is done by the DMaterialStepper class
576 
577 //}
578 
579 //---------------------------------
580 // GetCDCStereoWires
581 //---------------------------------
582 /// Extract the stereo wire data from the XML
583 bool DGeometry::GetCDCStereoWires(unsigned int ring,unsigned int ncopy,
584  double zcenter,double dz,
585  vector<vector<cdc_offset_t> >&cdc_offsets,
586  vector<DCDCWire*> &stereowires,
587  vector<double>&rot_angles,double dx,
588  double dy) const{
589  stringstream r_z_s,phi0_s,rot_s;
590 
591  // Create search strings for the straw geometrical properties
592  r_z_s << "//mposPhi[@volume='CDCstrawLong']/@R_Z/ring[@value='" << ring << "']";
593  phi0_s << "//mposPhi[@volume='CDCstrawLong']/@Phi0/ring[@value='" << ring << "']";
594  rot_s << "//mposPhi[@volume='CDCstrawLong']/@rot/ring[@value='" << ring << "']";
595 
596  vector<double>r_z;
597  double phi0;
598  vector<double>rot;
599 
600  // Extract data from the XML
601  if(!Get(r_z_s.str(), r_z)) return false;
602  if(!Get(phi0_s.str(), phi0)) return false;
603  if(!Get(rot_s.str(), rot)) return false;
604 
605  // Angular quantities
606  const double deg2rad=M_PI/180.;
607  double dphi=2*M_PI/double(ncopy);
608  phi0*=deg2rad;
609 
610  // Extract data from close-packed straws
611  double stereo=0.,stereo_sign=1.;
612  stereo=deg2rad*rot[0];
613  if (stereo<0.) stereo_sign=-1.;
614 
615  // Loop over the number of straws
616  for (unsigned int i=0;i<ncopy;i++){
617  DCDCWire *w=new DCDCWire;
618  double phi=phi0+double(i)*dphi;
619  w->ring=ring;
620  w->straw=i+1;
621 
622  // Find the nominal wire position and direction from the XML
623  DVector3 origin,udir;
624  origin.SetX(r_z[0]*cos(phi)+dx);
625  origin.SetY(r_z[0]*sin(phi)+dy);
626  origin.SetZ(zcenter);
627 
628  // Here, we need to define a coordinate system for the wire
629  // in which the wire runs along one axis. We call the directions
630  // of the axes in this coordinate system s,t, and u with
631  // the wire running in the "u" direction. The "s" direction
632  // will be defined by the direction pointing from the beamline
633  // to the midpoint of the wire.
634  udir.SetXYZ(0.0, 0.0,1.0);
635  udir.RotateX(stereo);
636  udir.RotateZ(phi);
637 
638  // Apply offsets in x and y
639  double half_dz=0.5*dz;
640  double x0=origin.x(),y0=origin.y();
641  double ux=udir.x()/udir.z();
642  double uy=udir.y()/udir.z();
643  unsigned int ringid=ring-1;
644  DVector3 downstream(x0+half_dz*ux+cdc_offsets[ringid][i].dx_d,
645  y0+half_dz*uy+cdc_offsets[ringid][i].dy_d,
646  zcenter+half_dz);
647  DVector3 upstream(x0-half_dz*ux+cdc_offsets[ringid][i].dx_u,
648  y0-half_dz*uy+cdc_offsets[ringid][i].dy_u,
649  zcenter-half_dz);
650  w->origin=0.5*(upstream+downstream);
651  w->origin.RotateX(rot_angles[0]);
652  w->origin.RotateY(rot_angles[1]);
653  w->origin.RotateZ(rot_angles[2]);
654 
655  w->phi=w->origin.Phi();
656 
657  // Wire direction
658  w->udir=downstream-upstream;
659  w->udir.RotateX(rot_angles[0]);
660  w->udir.RotateY(rot_angles[1]);
661  w->udir.RotateZ(rot_angles[2]);
662 
663  // For derivatives
664  w->udir_mag=w->udir.Mag();
665 
666  w->udir.SetMag(1.);
667  w->stereo=stereo_sign*w->udir.Angle(DVector3(0.,0.,1.));
668  // other directions for our wire coordinate system
669  w->sdir=w->origin;
670  w->sdir.SetMag(1.);
671  w->tdir = w->udir.Cross(w->sdir);
672 
673  // Some values needed for alignment derivatives
674  w->x0=dx; w->y0=dy; w->z0=zcenter;
675  w->phiX=rot_angles[0]; w->phiY=rot_angles[1]; w->phiZ=rot_angles[2];
676  w->r0=r_z[0]; w->phiStraw=phi; w->stereo_raw=stereo;
677 
678  stereowires.push_back(w);
679  }
680 
681  return true;
682 }
683 
684 //---------------------------------
685 // GetCDCAxialWires
686 //---------------------------------
687 /// Extract the axial wire data from the XML
688 bool DGeometry::GetCDCAxialWires(unsigned int ring,unsigned int ncopy,
689  double zcenter,double dz,
690  vector<vector<cdc_offset_t> >&cdc_offsets,
691  vector<DCDCWire*> &axialwires,
692  vector<double>&rot_angles,double dx,
693  double dy) const{
694  stringstream phi0_s,r_z_s;
695 
696  // Create search strings for the number of straws and the straw geometrical properties
697  phi0_s << "//mposPhi[@volume='CDCstrawShort']/@Phi0/ring[@value='" << ring << "']";
698  r_z_s << "//mposPhi[@volume='CDCstrawShort']/@R_Z/ring[@value='" << ring << "']";
699 
700  double phi0;
701  vector<double>r_z;
702 
703  // Extract the data from the XML
704  if(!Get(phi0_s.str(), phi0)) return false;
705  if(!Get(r_z_s.str(), r_z)) return false;
706 
707  // Angular quantities
708  double dphi=2*M_PI/double(ncopy);
709  phi0*=M_PI/180.;
710 
711  // Loop over the number of straws
712  for (unsigned int i=0;i<ncopy;i++){
713  DCDCWire *w=new DCDCWire;
714  double phi=phi0+double(i)*dphi;
715  w->ring=ring;
716  w->straw=i+1;
717 
718  // Find the nominal wire position from the XML
719  double x0=r_z[0]*cos(phi)+dx;
720  double y0=r_z[0]*sin(phi)+dy;
721 
722  // Apply offsets in x and y
723  double half_dz=0.5*dz;
724  unsigned int ringid=ring-1;
725  DVector3 downstream(x0+cdc_offsets[ringid][i].dx_d,
726  y0+cdc_offsets[ringid][i].dy_d,
727  zcenter+half_dz);
728  DVector3 upstream(x0+cdc_offsets[ringid][i].dx_u,
729  y0+cdc_offsets[ringid][i].dy_u,
730  zcenter-half_dz);
731  w->origin=0.5*(upstream+downstream);
732  w->origin.RotateX(rot_angles[0]);
733  w->origin.RotateY(rot_angles[1]);
734  w->origin.RotateZ(rot_angles[2]);
735  w->phi=w->origin.Phi();
736 
737  // Here, we need to define a coordinate system for the wire
738  // in which the wire runs along one axis. We call the directions
739  // of the axes in this coordinate system s,t, and u with
740  // the wire running in the "u" direction. The "s" direction
741  // will be defined by the direction pointing from the beamline
742  // to the midpoint of the wire.
743  w->udir=downstream-upstream;
744  w->udir.RotateX(rot_angles[0]);
745  w->udir.RotateY(rot_angles[1]);
746  w->udir.RotateZ(rot_angles[2]);
747 
748  // For derivatives
749  w->udir_mag=w->udir.Mag();
750 
751  w->udir.SetMag(1.);
752 
753  w->stereo=w->udir.Angle(DVector3(0.,0.,1.));
754  // other directions for our wire coordinate system
755  w->sdir=w->origin;
756  w->sdir.SetMag(1.);
757  w->tdir = w->udir.Cross(w->sdir);
758 
759  // Some values needed for alignment derivatives
760  w->x0=dx; w->y0=dy; w->z0=zcenter;
761  w->phiX=rot_angles[0]; w->phiY=rot_angles[1]; w->phiZ=rot_angles[2];
762  w->r0=r_z[0]; w->phiStraw=phi; w->stereo_raw=0.0;
763 
764  axialwires.push_back(w);
765  }
766 
767  return true;
768 }
769 
770 //---------------------------------
771 // GetCDCWires
772 //---------------------------------
773 bool DGeometry::GetCDCWires(vector<vector<DCDCWire *> >&cdcwires) const{
774  // Get nominal geometry from XML
775  vector<double>cdc_origin;
776  vector<double>cdc_length;
777  Get("//posXYZ[@volume='CentralDC']/@X_Y_Z",cdc_origin);
778  Get("//tubs[@name='STRW']/@Rio_Z",cdc_length);
779 
780  // Get CDC rotation angles
781  vector<double>rot_angles;
782  Get("//posXYZ[@volume='CentralDC']/@rot", rot_angles);
783  rot_angles[0]*=M_PI/180.;
784  rot_angles[1]*=M_PI/180.;
785  rot_angles[2]*=M_PI/180.;
786 
787  double dX=0.0, dY=0.0, dZ=0.0;
788  double dPhiX=0.0,dPhiY=0.0,dPhiZ=0.0;
789 
790  JCalibration * jcalib = dapp->GetJCalibration(runnumber);
791  vector<map<string,double> >vals;
792  if (jcalib->Get("CDC/global_alignment",vals)==false){
793  map<string,double> &row = vals[0];
794  dX=row["dX"];
795  dY=row["dY"];
796  dZ=row["dZ"];
797  dPhiX=row["dPhiX"];
798  dPhiY=row["dPhiY"];
799  dPhiZ=row["dPhiZ"];
800  }
801 
802 
803  // Shift the CDC origin according to alignment
804  cdc_origin[0]+=dX;
805  cdc_origin[1]+=dY;
806  cdc_origin[2]+=dZ;
807 
808  // Shift the CDC rotation according to the alignment
809  rot_angles[0]+=dPhiX;
810  rot_angles[1]+=dPhiY;
811  rot_angles[2]+=dPhiZ;
812 
813  double zmin=cdc_origin[2];
814  double zmax=zmin+cdc_length[2];
815  double zcenter=0.5*(zmin+zmax);
816  double L=zmax-zmin;
817 
818  // Get number of straws for each layer from the XML
819  unsigned int numstraws[28];
820  stringstream ncopy_s;
821 
822  // Get number of straws for each ring
823  for (unsigned int ring=1;ring<=28;ring++){
824  // Create search string for the number of straws
825  ncopy_s << "//section[@name='CentralDC']/composition/mposPhi/@ncopy/ring[@value='" << ring << "']";
826  Get(ncopy_s.str(),numstraws[ring-1]);
827  ncopy_s.str("");
828  ncopy_s.clear();
829  }
830 
831  // Get straw-by-straw offsets from calibration database
832  vector<cdc_offset_t>tempvec;
833  vector<vector<cdc_offset_t> >cdc_offsets;
834 
835  if (jcalib->Get("CDC/wire_alignment",vals)==false){
836  unsigned int straw_count=0,ring_count=0;
837  for(unsigned int i=0; i<vals.size(); i++){
838  map<string,double> &row = vals[i];
839 
840  // put the vector of offsets for the current ring into the offsets vector
841  if (straw_count==numstraws[ring_count]){
842  straw_count=0;
843  ring_count++;
844 
845  cdc_offsets.push_back(tempvec);
846 
847  tempvec.clear();
848  }
849 
850  // Get the offsets from the calibration database
852  temp.dx_u=row["dxu"];
853  //temp.dx_u=0.;
854 
855  temp.dy_u=row["dyu"];
856  //temp.dy_u=0.;
857 
858  temp.dx_d=row["dxd"];
859  //temp.dx_d=0.;
860 
861  temp.dy_d=row["dyd"];
862  //temp.dy_d=0.;
863 
864  tempvec.push_back(temp);
865 
866  straw_count++;
867  }
868  cdc_offsets.push_back(tempvec);
869  }
870  else{
871  jerr<< "CDC wire alignment table not available... bailing... " <<endl;
872  exit(0);
873  }
874 
875  // First axial layer
876  for (unsigned int ring=1;ring<5;ring++){
877  vector<DCDCWire*>straws;
878  if (!GetCDCAxialWires(ring,numstraws[ring-1],zcenter,L,cdc_offsets,straws,
879  rot_angles,cdc_origin[0],cdc_origin[1])) return false;
880  cdcwires.push_back(straws);
881  }
882 
883  // First set of stereo layers
884  for (unsigned int i=0;i<8;i++){
885  vector<DCDCWire*>straws;
886  if (!GetCDCStereoWires(i+5,numstraws[i+4],zcenter,L,cdc_offsets,straws,
887  rot_angles,cdc_origin[0],cdc_origin[1])) return false;
888  cdcwires.push_back(straws);
889  }
890 
891  // Second axial layer
892  for (unsigned int ring=13;ring<17;ring++){
893  vector<DCDCWire*>straws;
894  if (!GetCDCAxialWires(ring,numstraws[ring-1],zcenter,L,cdc_offsets,straws,
895  rot_angles,cdc_origin[0],cdc_origin[1])) return false;
896  cdcwires.push_back(straws);
897  }
898 
899  // Second set of stereo layers
900  for (unsigned int i=8;i<16;i++){
901  vector<DCDCWire*>straws;
902  if (!GetCDCStereoWires(i+9,numstraws[i+8],zcenter,L,cdc_offsets,straws,
903  rot_angles,cdc_origin[0],cdc_origin[1])) return false;
904  cdcwires.push_back(straws);
905  }
906 
907  // Third axial layer
908  for (unsigned int ring=25;ring<29;ring++){
909  vector<DCDCWire*>straws;
910  if (!GetCDCAxialWires(ring,numstraws[ring-1],zcenter,L,cdc_offsets,straws,
911  rot_angles,cdc_origin[0],cdc_origin[1])) return false;
912  cdcwires.push_back(straws);
913  }
914 
915  // Calculate wire lengths and compute "s" and "t" direction vectors (orthogonal to "u")
916  for (unsigned int i=0;i<cdcwires.size();i++){
917  for (unsigned int j=0;j<cdcwires[i].size();j++){
918  DCDCWire *w=cdcwires[i][j];
919  w->L=L/cos(w->stereo);
920 
921  // With the addition of close-packed stereo wires, the vector connecting
922  // the center of the wire to the beamline ("s" direction) is not necessarily
923  // perpendicular to the beamline. By definition, we want the "s" direction
924  // to be perpendicular to the wire direction "u" and pointing at the beamline.
925  //
926  // NOTE: This extensive comment is here because the result, when implmented
927  // below caused a WORSE residual distribution in the close-packed stereo
928  // layers. Owing to lack of time currently to track the issue down (most
929  // likely in DReferenceTrajectory) I'm commenting out the "correct" calculation
930  // of s, but leaving this comment so the issue can be revisited later. This
931  // error leads to around 100 micron errors in the C.P.S. wires, but they
932  // are completely washed out when the position smearing of 150 microns is applied
933  // making the error unnoticable except when position smearing is not applied.
934  //
935  // April 2, 2009 D.L.
936  //
937  // Here is how this is calculated -- We define a vector equation with 2 unknowns
938  // Z and S:
939  //
940  // Zz + Ss = W
941  //
942  // where: z = unit vector in z direction
943  // s = unit vector in "s" direction
944  // W = vector pointing to center of wire in lab coordinates
945  //
946  // Rearranging, we get:
947  //
948  // s = (W - Zz)/S
949  //
950  // Since s must be perpendicular to u, we take a dot product of s and u
951  // and set it equal to zero to determine Z:
952  //
953  // u.s = 0 = u.(W - Zz)/S => u.W = Zu.z
954  //
955  // or
956  //
957  // Z = u.W/u.z
958  //
959  // Thus, the s direction is just given by (W - (u.W/u.z)z)
960  //
961 
962  //w->sdir=w->origin-DVector3(0,0,w->origin.Z());
963  w->sdir = w->origin - DVector3(0.0, 0.0, w->udir.Dot(w->origin)/w->udir.Z()); // see above comments
964  w->sdir.SetMag(1.0);
965 
966  w->tdir = w->udir.Cross(w->sdir);
967  w->tdir.SetMag(1.0); // This isn't really needed
968  }
969  }
970 
971  return true;
972 }
973 
974 
975 //---------------------------------
976 // GetFDCCathodes
977 //---------------------------------
978 bool DGeometry::GetFDCCathodes(vector<vector<DFDCCathode *> >&fdccathodes) const{
979  // Get offsets tweaking nominal geometry from calibration database
980  JCalibration * jcalib = dapp->GetJCalibration(runnumber);
981  vector<map<string,double> >vals;
982  vector<fdc_cathode_offset_t>fdc_cathode_offsets;
983  if (jcalib->Get("FDC/cathode_alignment",vals)==false){
984  for(unsigned int i=0; i<vals.size(); i++){
985  map<string,double> &row = vals[i];
986 
987  // Get the offsets from the calibration database
989  temp.du=row["dU"];
990  //temp.du=0.;
991 
992  temp.dphi=row["dPhiU"];
993  //temp.dphi=0.;
994 
995  fdc_cathode_offsets.push_back(temp);
996 
997  temp.du=row["dV"];
998  //temp.du=0.;
999 
1000  temp.dphi=row["dPhiV"];
1001  //temp.dphi=0.;
1002 
1003  fdc_cathode_offsets.push_back(temp);
1004  }
1005  }
1006  vector< vector<double> >fdc_cathode_pitches;
1007  if (jcalib->Get("FDC/strip_pitches_v2",vals)==false){
1008  for(unsigned int i=0; i<vals.size(); i++){
1009  map<string,double> &row = vals[i];
1010 
1011  vector<double> uvals;
1012  // Get the offsets from the calibration database
1013  uvals.push_back(row["U_SP_1"]);
1014  uvals.push_back(row["U_G_1"]);
1015  uvals.push_back(row["U_SP_2"]);
1016  uvals.push_back(row["U_G_2"]);
1017  uvals.push_back(row["U_SP_3"]);
1018 
1019  fdc_cathode_pitches.push_back(uvals);
1020 
1021  vector<double> vvals;
1022  // Get the offsets from the calibration database
1023  vvals.push_back(row["V_SP_1"]);
1024  vvals.push_back(row["V_G_1"]);
1025  vvals.push_back(row["V_SP_2"]);
1026  vvals.push_back(row["V_G_2"]);
1027  vvals.push_back(row["V_SP_3"]);
1028 
1029  fdc_cathode_pitches.push_back(vvals);
1030  }
1031  }
1032  else{
1033  jerr << "Strip pitch calibration unavailable -- setting default..." <<endl;
1034  // set some sensible default
1035  for (unsigned int i=0;i<2*FDC_NUM_LAYERS;i++){
1036 
1037  vector<double> val;
1038  for (int j = 0; j < 5; j++){
1039  val.push_back(STRIP_SPACING);
1040  }
1041 
1042  fdc_cathode_pitches.push_back(val);
1043  }
1044  }
1045 
1046  // Generate the vector of cathode plane parameters
1047  for (int i=0;i<2*FDC_NUM_LAYERS; i++){
1048  double angle=(i%2)?(M_PI-CATHODE_ROT_ANGLE):(CATHODE_ROT_ANGLE);
1049 
1050  angle+=fdc_cathode_offsets[i].dphi;
1051  double SP1 = fdc_cathode_pitches[i][0];
1052  double SG1 = fdc_cathode_pitches[i][1];
1053  double SP2 = fdc_cathode_pitches[i][2];
1054  double SG2 = fdc_cathode_pitches[i][3];
1055  double SP3 = fdc_cathode_pitches[i][4];
1056 
1057  vector<DFDCCathode *>temp;
1058  for (int j=0; j<STRIPS_PER_PLANE; j++){
1059  DFDCCathode *c = new DFDCCathode;
1060  c->layer = i+1;
1061  c->strip = j+1;
1062  c->angle = angle;
1063  if (j<48) c->u=(-47.5*SP2 - SG1 + (j-47)*SP1) + fdc_cathode_offsets[i].du;
1064  else if (j<144) c->u=(double(j)-95.5)*SP2 + fdc_cathode_offsets[i].du;
1065  else c->u=(47.5*SP2 + SG2 + (j-144)*SP3) + fdc_cathode_offsets[i].du;
1066 
1067  temp.push_back(c);
1068  }
1069  fdccathodes.push_back(temp);
1070  }
1071 
1072  return true;
1073 }
1074 
1075 //---------------------------------
1076 // GetFDCWires
1077 //---------------------------------
1078 bool DGeometry::GetFDCWires(vector<vector<DFDCWire *> >&fdcwires) const{
1079  // Get geometrical information from database
1080  vector<double>z_wires;
1081  vector<double>stereo_angles;
1082 
1083  if(!GetFDCZ(z_wires)) return false;
1084  if(!GetFDCStereo(stereo_angles)) return false;
1085 
1086  // Get package rotation angles
1087  double ThetaX[4],ThetaY[4],ThetaZ[4];
1088  vector<double>rot_angles;
1089  Get("//posXYZ[@volume='forwardDC_package_1']/@rot", rot_angles);
1090  ThetaX[0]=rot_angles[0]*M_PI/180.;
1091  ThetaY[0]=rot_angles[1]*M_PI/180.;
1092  ThetaZ[0]=rot_angles[2]*M_PI/180.;
1093  Get("//posXYZ[@volume='forwardDC_package_2']/@rot", rot_angles);
1094  ThetaX[1]=rot_angles[0]*M_PI/180.;
1095  ThetaY[1]=rot_angles[1]*M_PI/180.;
1096  ThetaZ[1]=rot_angles[2]*M_PI/180.;
1097  Get("//posXYZ[@volume='forwardDC_package_3']/@rot", rot_angles);
1098  ThetaX[2]=rot_angles[0]*M_PI/180.;
1099  ThetaY[2]=rot_angles[1]*M_PI/180.;
1100  ThetaZ[2]=rot_angles[2]*M_PI/180.;
1101  Get("//posXYZ[@volume='forwardDC_package_4']/@rot", rot_angles);
1102  ThetaX[3]=rot_angles[0]*M_PI/180.;
1103  ThetaY[3]=rot_angles[1]*M_PI/180.;
1104  ThetaZ[3]=rot_angles[2]*M_PI/180.;
1105  // Get package offsets
1106  double dX[4],dY[4];
1107  vector<double>offsets;
1108  Get("//posXYZ[@volume='forwardDC_package_1']/@X_Y_Z",offsets);
1109  dX[0]=offsets[0];
1110  dY[0]=offsets[1];
1111  Get("//posXYZ[@volume='forwardDC_package_2']/@X_Y_Z",offsets);
1112  dX[1]=offsets[0];
1113  dY[1]=offsets[1];
1114  Get("//posXYZ[@volume='forwardDC_package_3']/@X_Y_Z",offsets);
1115  dX[2]=offsets[0];
1116  dY[2]=offsets[1];
1117  Get("//posXYZ[@volume='forwardDC_package_4']/@X_Y_Z",offsets);
1118  dX[3]=offsets[0];
1119  dY[3]=offsets[1];
1120 
1121  // Get offsets tweaking nominal geometry from calibration database
1122  JCalibration * jcalib = dapp->GetJCalibration(runnumber);
1123  vector<map<string,double> >vals;
1124  vector<fdc_wire_offset_t>fdc_wire_offsets;
1125  if (jcalib->Get("FDC/wire_alignment",vals)==false){
1126  for(unsigned int i=0; i<vals.size(); i++){
1127  map<string,double> &row = vals[i];
1128 
1129  // Get the offsets from the calibration database
1131  temp.du=row["dU"];
1132  //temp.du=0.;
1133 
1134  temp.dphi=row["dPhi"];
1135  //temp.dphi=0.;
1136 
1137  temp.dz=row["dZ"];
1138  // temp.dz=0.;
1139 
1140  fdc_wire_offsets.push_back(temp);
1141  }
1142  }
1143 
1144  vector<fdc_wire_rotation_t>fdc_wire_rotations;
1145  if (jcalib->Get("FDC/cell_rotations",vals)==false){
1146  for(unsigned int i=0; i<vals.size(); i++){
1147  map<string,double> &row = vals[i];
1148 
1149  // Get the offsets from the calibration database
1151  temp.dPhiX=row["dPhiX"];
1152  temp.dPhiY=row["dPhiY"];
1153  temp.dPhiZ=row["dPhiZ"];
1154 
1155  fdc_wire_rotations.push_back(temp);
1156  }
1157  }
1158 
1159  // Generate the vector of wire plane parameters
1160  for(int i=0; i<FDC_NUM_LAYERS; i++){
1161  double angle=-stereo_angles[i]*M_PI/180.+fdc_wire_offsets[i].dphi;
1162 
1163  vector<DFDCWire *>temp;
1164  for(int j=0; j<WIRES_PER_PLANE; j++){
1165  unsigned int pack_id=i/6;
1166 
1167  DFDCWire *w = new DFDCWire;
1168  w->layer = i+1;
1169  w->wire = j+1;
1170  w->angle = angle;
1171 
1172  // find coordinates of center of wire in rotated system
1173  float u = U_OF_WIRE_ZERO + WIRE_SPACING*(float)(j);
1174  w->u=u+fdc_wire_offsets[i].du;
1175 
1176  // Rotate coordinates into lab system and set the wire's origin
1177  // Note that the FDC measures "angle" such that angle=0
1178  // corresponds to the anode wire in the vertical direction
1179  // (i.e. at phi=90 degrees).
1180  float x = u*sin(angle + M_PI/2.0);
1181  float y = u*cos(angle + M_PI/2.0);
1182  w->origin.SetXYZ(x,y,0.);
1183  w->origin.RotateX(ThetaX[pack_id]+fdc_wire_rotations[i].dPhiX);
1184  w->origin.RotateY(ThetaY[pack_id]+fdc_wire_rotations[i].dPhiY);
1185  w->origin.RotateZ(ThetaZ[pack_id]+fdc_wire_rotations[i].dPhiZ);
1186  DVector3 globalOffsets(dX[pack_id],dY[pack_id],z_wires[i]+fdc_wire_offsets[i].dz);
1187  w->origin+=globalOffsets;
1188 
1189  // Length of wire is set by active radius
1190  w->L = 2.0*sqrt(pow(FDC_ACTIVE_RADIUS,2.0) - u*u);
1191 
1192  // Set directions of wire's coordinate system with "udir"
1193  // along wire.
1194  w->udir.SetXYZ(sin(angle),cos(angle),0.0);
1195  w->udir.RotateX(ThetaX[pack_id]+fdc_wire_rotations[i].dPhiX);
1196  w->udir.RotateY(ThetaY[pack_id]+fdc_wire_rotations[i].dPhiY);
1197  w->udir.RotateZ(ThetaZ[pack_id]+fdc_wire_rotations[i].dPhiZ);
1198  w->angles.SetXYZ(ThetaX[pack_id]+fdc_wire_rotations[i].dPhiX,
1199  ThetaY[pack_id]+fdc_wire_rotations[i].dPhiY,
1200  ThetaZ[pack_id]+fdc_wire_rotations[i].dPhiZ);
1201  w->u+=dX[pack_id]*w->udir.y()-dY[pack_id]*w->udir.x();
1202 
1203  // "s" points in direction from beamline to midpoint of
1204  // wire. This happens to be the same direction as "origin"
1205  w->sdir = w->origin;
1206  w->sdir.SetMag(1.0);
1207 
1208  w->tdir = w->udir.Cross(w->sdir);
1209  w->tdir.SetMag(1.0); // This isn't really needed
1210  temp.push_back(w);
1211  }
1212  fdcwires.push_back(temp);
1213  }
1214 
1215  return true;
1216 }
1217 
1218 //---------------------------------
1219 // GetFDCZ
1220 //---------------------------------
1221 bool DGeometry::GetFDCZ(vector<double> &z_wires) const
1222 {
1223  // The FDC geometry is defined as 4 packages, each containing 2
1224  // "module"s and each of those containing 3 "chambers". The modules
1225  // are placed as multiple copies in Z using mposZ, but none of the
1226  // others are (???).
1227  //
1228  // This method is currently hardwired to assume 4 packages and
1229  // 3 chambers. (The number of modules is discovered via the
1230  // "ncopy" attribute of mposZ.)
1231 
1232  vector<double> ForwardDC;
1233  vector<double> forwardDC;
1234  vector<double> forwardDC_package[4];
1235  vector<double> forwardDC_module[4];
1236  vector<double> forwardDC_chamber[4][6];
1237 
1238  if(!Get("//section/composition/posXYZ[@volume='ForwardDC']/@X_Y_Z", ForwardDC)) return false;
1239  if(!Get("//composition[@name='ForwardDC']/posXYZ[@volume='forwardDC']/@X_Y_Z", forwardDC)) return false;
1240  if(!Get("//posXYZ[@volume='forwardDC_package_1']/@X_Y_Z", forwardDC_package[0])) return false;
1241  if(!Get("//posXYZ[@volume='forwardDC_package_2']/@X_Y_Z", forwardDC_package[1])) return false;
1242  if(!Get("//posXYZ[@volume='forwardDC_package_3']/@X_Y_Z", forwardDC_package[2])) return false;
1243  if(!Get("//posXYZ[@volume='forwardDC_package_4']/@X_Y_Z", forwardDC_package[3])) return false;
1244  if(!Get("//posXYZ[@volume='forwardDC_module_1']/@X_Y_Z", forwardDC_module[0])) return false;
1245  if(!Get("//posXYZ[@volume='forwardDC_module_2']/@X_Y_Z", forwardDC_module[1])) return false;
1246  if(!Get("//posXYZ[@volume='forwardDC_module_3']/@X_Y_Z", forwardDC_module[2])) return false;
1247  if(!Get("//posXYZ[@volume='forwardDC_module_4']/@X_Y_Z", forwardDC_module[3])) return false;
1248  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='1']", forwardDC_chamber[0][0])) return false;
1249  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='2']", forwardDC_chamber[0][1])) return false;
1250  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='3']", forwardDC_chamber[0][2])) return false;
1251  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='4']", forwardDC_chamber[0][3])) return false;
1252  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='5']", forwardDC_chamber[0][4])) return false;
1253  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@X_Y_Z/layer[@value='6']", forwardDC_chamber[0][5])) return false;
1254  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='1']", forwardDC_chamber[1][0])) return false;
1255  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='2']", forwardDC_chamber[1][1])) return false;
1256  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='3']", forwardDC_chamber[1][2])) return false;
1257  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='4']", forwardDC_chamber[1][3])) return false;
1258  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='5']", forwardDC_chamber[1][4])) return false;
1259  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@X_Y_Z/layer[@value='6']", forwardDC_chamber[1][5])) return false;
1260  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='1']", forwardDC_chamber[2][0])) return false;
1261  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='2']", forwardDC_chamber[2][1])) return false;
1262  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='3']", forwardDC_chamber[2][2])) return false;
1263  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='4']", forwardDC_chamber[2][3])) return false;
1264  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='5']", forwardDC_chamber[2][4])) return false;
1265  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@X_Y_Z/layer[@value='6']", forwardDC_chamber[2][5])) return false;
1266  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='1']", forwardDC_chamber[3][0])) return false;
1267  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='2']", forwardDC_chamber[3][1])) return false;
1268  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='3']", forwardDC_chamber[3][2])) return false;
1269  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='4']", forwardDC_chamber[3][3])) return false;
1270  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='5']", forwardDC_chamber[3][4])) return false;
1271  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@X_Y_Z/layer[@value='6']", forwardDC_chamber[3][5])) return false;
1272 
1273  // Offset due to global FDC envelopes
1274  double zfdc = ForwardDC[2] + forwardDC[2];
1275 
1276  // Loop over packages
1277  for(int package=1; package<=4; package++){
1278  double z_package = forwardDC_package[package-1][2];
1279 
1280  // Each "package" has 1 "module" which is currently positioned at 0,0,0 but
1281  // that could in principle change so we add the module z-offset
1282  double z_module = forwardDC_module[package-1][2];
1283 
1284  // Loop over chambers in this module
1285  for(int chamber=1; chamber<=6; chamber++){
1286  double z_chamber = forwardDC_chamber[package-1][chamber-1][2];
1287 
1288  double z = zfdc + z_package + z_module + z_chamber;
1289  z_wires.push_back(z);
1290  }
1291  }
1292 
1293  return true;
1294 }
1295 
1296 //---------------------------------
1297 // GetFDCStereo
1298 //---------------------------------
1299 bool DGeometry::GetFDCStereo(vector<double> &stereo_angles) const
1300 {
1301  // The FDC geometry is defined as 4 packages, each containing 2
1302  // "module"s and each of those containing 3 "chambers". The modules
1303  // are placed as multiple copies in Z using mposZ, but none of the
1304  // others are (???).
1305  //
1306  // This method is currently hardwired to assume 4 packages and
1307  // 3 chambers. (The number of modules is discovered via the
1308  // "ncopy" attribute of mposZ.)
1309  //
1310  // Stereo angles are assumed to be rotated purely about the z-axis
1311  // and the units are not specified, but the XML currently uses degrees.
1312 
1313  vector<double> forwardDC_module[4];
1314  vector<double> forwardDC_chamber[4][6];
1315 
1316  if(!Get("//posXYZ[@volume='forwardDC_module_1']/@X_Y_Z", forwardDC_module[0])) return false;
1317  if(!Get("//posXYZ[@volume='forwardDC_module_2']/@X_Y_Z", forwardDC_module[1])) return false;
1318  if(!Get("//posXYZ[@volume='forwardDC_module_3']/@X_Y_Z", forwardDC_module[2])) return false;
1319  if(!Get("//posXYZ[@volume='forwardDC_module_4']/@X_Y_Z", forwardDC_module[3])) return false;
1320  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='1']", forwardDC_chamber[0][0])) return false;
1321  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='2']", forwardDC_chamber[0][1])) return false;
1322  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='3']", forwardDC_chamber[0][2])) return false;
1323  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='4']", forwardDC_chamber[0][3])) return false;
1324  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='5']", forwardDC_chamber[0][4])) return false;
1325  if(!Get("//posXYZ[@volume='forwardDC_chamber_1']/@rot/layer[@value='6']", forwardDC_chamber[0][5])) return false;
1326  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='1']", forwardDC_chamber[1][0])) return false;
1327  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='2']", forwardDC_chamber[1][1])) return false;
1328  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='3']", forwardDC_chamber[1][2])) return false;
1329  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='4']", forwardDC_chamber[1][3])) return false;
1330  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='5']", forwardDC_chamber[1][4])) return false;
1331  if(!Get("//posXYZ[@volume='forwardDC_chamber_2']/@rot/layer[@value='6']", forwardDC_chamber[1][5])) return false;
1332  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='1']", forwardDC_chamber[2][0])) return false;
1333  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='2']", forwardDC_chamber[2][1])) return false;
1334  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='3']", forwardDC_chamber[2][2])) return false;
1335  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='4']", forwardDC_chamber[2][3])) return false;
1336  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='5']", forwardDC_chamber[2][4])) return false;
1337  if(!Get("//posXYZ[@volume='forwardDC_chamber_3']/@rot/layer[@value='6']", forwardDC_chamber[2][5])) return false;
1338  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='1']", forwardDC_chamber[3][0])) return false;
1339  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='2']", forwardDC_chamber[3][1])) return false;
1340  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='3']", forwardDC_chamber[3][2])) return false;
1341  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='4']", forwardDC_chamber[3][3])) return false;
1342  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='5']", forwardDC_chamber[3][4])) return false;
1343  if(!Get("//posXYZ[@volume='forwardDC_chamber_4']/@rot/layer[@value='6']", forwardDC_chamber[3][5])) return false;
1344 
1345  // Loop over packages
1346  for(int package=1; package<=4; package++){
1347 
1348  // Loop over chambers
1349  for(int chamber=1; chamber<=6; chamber++){
1350  // if (chamber==4) forwardDC_chamber[package-1][chamber-1][2]+=15.0;
1351  stereo_angles.push_back(forwardDC_chamber[package-1][chamber-1][2]);
1352  }
1353  }
1354 
1355  return true;
1356 }
1357 
1358 //---------------------------------
1359 // GetFDCRmin
1360 //---------------------------------
1361 bool DGeometry::GetFDCRmin(vector<double> &rmin_packages) const
1362 {
1363  vector<double> FDA[4];
1364 
1365  if(!Get("//section[@name='ForwardDC']/tubs[@name='FDA1']/@Rio_Z", FDA[0])) return false;
1366  if(!Get("//section[@name='ForwardDC']/tubs[@name='FDA2']/@Rio_Z", FDA[1])) return false;
1367  if(!Get("//section[@name='ForwardDC']/tubs[@name='FDA3']/@Rio_Z", FDA[2])) return false;
1368  if(!Get("//section[@name='ForwardDC']/tubs[@name='FDA4']/@Rio_Z", FDA[3])) return false;
1369 
1370  rmin_packages.push_back(FDA[0][0]);
1371  rmin_packages.push_back(FDA[1][0]);
1372  rmin_packages.push_back(FDA[2][0]);
1373  rmin_packages.push_back(FDA[3][0]);
1374 
1375  return true;
1376 }
1377 
1378 //---------------------------------
1379 // GetFDCRmax
1380 //---------------------------------
1381 bool DGeometry::GetFDCRmax(double &rmax_active_fdc) const
1382 {
1383  // We assume that all packages have the same outer radius of the
1384  // active area.
1385  vector<double> FDA1;
1386 
1387  bool good = Get("//section[@name='ForwardDC']/tubs[@name='FDA1']/@Rio_Z", FDA1);
1388 
1389  if(!good){
1390  _DBG_<<"Unable to retrieve FDC Rmax values."<<endl;
1391  return good;
1392  }
1393 
1394  rmax_active_fdc = FDA1[1];
1395 
1396  return good;
1397 }
1398 
1399 //---------------------------------
1400 // GetCDCOption
1401 //---------------------------------
1402 bool DGeometry::GetCDCOption(string &cdc_option) const
1403 {
1404  bool good = Get("//section[@name='CentralDC_s']/composition/posXYZ/@volume", cdc_option);
1405 
1406  if(!good){
1407  _DBG_<<"Unable to retrieve CDC option string."<<endl;
1408  }
1409 
1410  return good;
1411 }
1412 
1413 //---------------------------------
1414 // GetCDCCenterZ
1415 //---------------------------------
1416 bool DGeometry::GetCDCCenterZ(double &cdc_center_z) const
1417 {
1418 
1419  return false;
1420 }
1421 
1422 //---------------------------------
1423 // GetCDCAxialLength
1424 //---------------------------------
1425 bool DGeometry::GetCDCAxialLength(double &cdc_axial_length) const
1426 {
1427  vector<double> Rio_Z;
1428  bool good = Get("//section[@name='CentralDC']/tubs[@name='STRW']/@Rio_Z", Rio_Z);
1429  cdc_axial_length = Rio_Z[2];
1430 
1431  if(!good){
1432  _DBG_<<"Unable to retrieve CDC axial wire length"<<endl;
1433  }
1434 
1435  return false;
1436 }
1437 
1438 //---------------------------------
1439 // GetCDCStereo
1440 //---------------------------------
1441 bool DGeometry::GetCDCStereo(vector<double> &cdc_stereo) const
1442 {
1443 
1444  return false;
1445 }
1446 
1447 //---------------------------------
1448 // GetCDCRmid
1449 //---------------------------------
1450 bool DGeometry::GetCDCRmid(vector<double> &cdc_rmid) const
1451 {
1452 
1453  return false;
1454 }
1455 
1456 //---------------------------------
1457 // GetCDCNwires
1458 //---------------------------------
1459 bool DGeometry::GetCDCNwires(vector<int> &cdc_nwires) const
1460 {
1461  cdc_nwires.push_back(42);
1462  cdc_nwires.push_back(42);
1463  cdc_nwires.push_back(54);
1464  cdc_nwires.push_back(54);
1465  cdc_nwires.push_back(66);
1466  cdc_nwires.push_back(66);
1467  cdc_nwires.push_back(80);
1468  cdc_nwires.push_back(80);
1469  cdc_nwires.push_back(93);
1470  cdc_nwires.push_back(93);
1471  cdc_nwires.push_back(106);
1472  cdc_nwires.push_back(106);
1473  cdc_nwires.push_back(123);
1474  cdc_nwires.push_back(123);
1475  cdc_nwires.push_back(135);
1476  cdc_nwires.push_back(135);
1477  cdc_nwires.push_back(146);
1478  cdc_nwires.push_back(146);
1479  cdc_nwires.push_back(158);
1480  cdc_nwires.push_back(158);
1481  cdc_nwires.push_back(170);
1482  cdc_nwires.push_back(170);
1483  cdc_nwires.push_back(182);
1484  cdc_nwires.push_back(182);
1485  cdc_nwires.push_back(197);
1486  cdc_nwires.push_back(197);
1487  cdc_nwires.push_back(209);
1488  cdc_nwires.push_back(209);
1489 
1490  return false;
1491 }
1492 
1493 
1494 //---------------------------------
1495 // GetCDCEndplate
1496 //---------------------------------
1497 bool DGeometry::GetCDCEndplate(double &z,double &dz,double &rmin,double &rmax)
1498  const{
1499 
1500  vector<double>cdc_origin;
1501  vector<double>cdc_center;
1502  vector<double>cdc_layers_offset;
1503  vector<double>cdc_endplate_pos;
1504  vector<double>cdc_endplate_dim;
1505 
1506  if(!Get("//posXYZ[@volume='CentralDC'/@X_Y_Z",cdc_origin)) return false;
1507  if(!Get("//posXYZ[@volume='centralDC']/@X_Y_Z",cdc_center)) return false;
1508  if(!Get("//posXYZ[@volume='CDPD']/@X_Y_Z",cdc_endplate_pos)) return false;
1509  if(!Get("//tubs[@name='CDPD']/@Rio_Z",cdc_endplate_dim)) return false;
1510  if(!Get("//posXYZ[@volume='CDClayers']/@X_Y_Z",cdc_layers_offset)) return false;
1511 
1512  if(cdc_origin.size()<3){
1513  _DBG_<<"cdc_origin.size()<3 !"<<endl;
1514  return false;
1515  }
1516  if(cdc_center.size()<3){
1517  _DBG_<<"cdc_center.size()<3 !"<<endl;
1518  return false;
1519  }
1520  if(cdc_endplate_pos.size()<3){
1521  _DBG_<<"cdc_endplate_pos.size()<3 !"<<endl;
1522  return false;
1523  }
1524  if(cdc_endplate_dim.size()<3){
1525  _DBG_<<"cdc_endplate_dim.size()<3 !"<<endl;
1526  return false;
1527  }
1528  if (cdc_layers_offset.size()<3){
1529  _DBG_<<"cdc_layers_offset.size()<3 !"<<endl;
1530  return false;
1531  }
1532 
1533  z=cdc_origin[2]+cdc_center[2]+cdc_endplate_pos[2]+cdc_layers_offset[2];
1534  dz=cdc_endplate_dim[2];
1535  rmin=cdc_endplate_dim[0];
1536  rmax=cdc_endplate_dim[1];
1537 
1538  return true;
1539  }
1540 //---------------------------------
1541 // GetBCALRmin
1542 //---------------------------------
1543 // Including the support plate
1544 bool DGeometry::GetBCALRmin(float &bcal_rmin) const
1545 {
1546  vector<float> bcal_mother_Rio_Z;
1547  bool good = Get("//section[@name='BarrelEMcal']/tubs[@name='BCAL']/@Rio_Z", bcal_mother_Rio_Z);
1548  if(!good){
1549  _DBG_<<"Unable to retrieve BCAL mother RioZ info."<<endl;
1550  bcal_rmin = 0.0;
1551  return false;
1552  }
1553  if(bcal_mother_Rio_Z.size() == 3){
1554  bcal_rmin = bcal_mother_Rio_Z[0];
1555  return true;
1556  }
1557  else{
1558  _DBG_<<"Wrong vector size for BCAL mother RioZ!!!"<<endl;
1559  bcal_rmin = 0.0;
1560  return false;
1561  }
1562 }
1563 
1564 //---------------------------------
1565 // GetBCALfADCRadii
1566 //---------------------------------
1567 bool DGeometry::GetBCALfADCRadii(vector<float> &fADC_radii) const
1568 {
1569  vector<float> BM[5];
1570 
1571  if(!Get("//section[@name='BarrelEMcal']/tubs[@name='BM01']/@Rio_Z", BM[0])) return false;
1572  if(!Get("//section[@name='BarrelEMcal']/tubs[@name='BM02']/@Rio_Z", BM[1])) return false;
1573  if(!Get("//section[@name='BarrelEMcal']/tubs[@name='BM04']/@Rio_Z", BM[2])) return false;
1574  if(!Get("//section[@name='BarrelEMcal']/tubs[@name='BMF7']/@Rio_Z", BM[3])) return false;
1575  if(!Get("//section[@name='BarrelEMcal']/tubs[@name='BMFA']/@Rio_Z", BM[4])) return false;
1576 
1577  fADC_radii.push_back(BM[0][0]);
1578  fADC_radii.push_back(BM[1][0]);
1579  fADC_radii.push_back(BM[2][0]);
1580  fADC_radii.push_back(BM[3][0]);
1581  fADC_radii.push_back(BM[4][1]);
1582 
1583  return true;
1584 }
1585 
1586 //---------------------------------
1587 // GetBCALNmodules
1588 //---------------------------------
1589 bool DGeometry::GetBCALNmodules(unsigned int &bcal_nmodules) const
1590 {
1591  vector<unsigned int> ncopy;
1592  bool good = Get("//section[@name='BarrelEMcal']/composition/mposPhi/@ncopy", ncopy);
1593  if(!good){
1594  _DBG_<<"Unable to retrieve BCAL barrelModule ncopy info."<<endl;
1595  bcal_nmodules = 0;
1596  return false;
1597  }
1598  if(ncopy.size() == 1){
1599  bcal_nmodules = ncopy[0];
1600  return true;
1601  }else{
1602  _DBG_<<"Wrong vector size for BCAL barrelModule ncopy!!!"<<endl;
1603  bcal_nmodules = 0;
1604  return false;
1605  }
1606 }
1607 
1608 //---------------------------------
1609 // GetBCALCenterZ
1610 //---------------------------------
1611 bool DGeometry::GetBCALCenterZ(float &bcal_center_z) const
1612 {
1613  vector<float> z0;
1614  bool good = Get("//section[@name='BarrelEMcal']/parameters/real[@name='z0']/@value", z0);
1615  if(!good){
1616  _DBG_<<"Unable to retrieve BCAL parameters z0 info."<<endl;
1617  bcal_center_z = 0.0;
1618  return false;
1619  }
1620  if(z0.size() == 1){
1621  bcal_center_z = z0[0];
1622  return true;
1623  }else{
1624  _DBG_<<"Wrong vector size for BCAL parameters z0!!!"<<endl;
1625  bcal_center_z = 0.0;
1626  return false;
1627  }
1628 
1629 }
1630 
1631 //---------------------------------
1632 // GetBCALLength
1633 //---------------------------------
1634 // The lightguides are not included
1635 bool DGeometry::GetBCALLength(float &bcal_length) const
1636 {
1637  vector<float> module_length;
1638  bool good = Get("//section[@name='BarrelEMcal']/tubs[@name='BM01']/@Rio_Z", module_length);
1639  if(!good){
1640  _DBG_<<"Unable to retrieve BCAL submodule RioZ info."<<endl;
1641  bcal_length = 0.0;
1642  return false;
1643  }
1644  if(module_length.size() == 3){
1645  bcal_length = module_length[2];
1646  return true;
1647  }
1648  else{
1649  _DBG_<<"Wrong vector size for BCAL submodule RioZ!!!"<<endl;
1650  bcal_length = 0.0;
1651  return false;
1652  }
1653 }
1654 
1655 //---------------------------------
1656 // GetBCALDepth
1657 //---------------------------------
1658 // Including the support plate and the support bar
1659 bool DGeometry::GetBCALDepth(float &bcal_depth) const
1660 {
1661  vector<float> bcal_moth_Rio_Z;
1662  bool good = Get("//section[@name='BarrelEMcal']/tubs[@name='BCAL']/@Rio_Z", bcal_moth_Rio_Z);
1663  if(!good){
1664  _DBG_<<"Unable to retrieve BCAL mother RioZ info."<<endl;
1665  bcal_depth = 0.0;
1666  return false;
1667  }
1668  if(bcal_moth_Rio_Z.size() == 3){
1669  bcal_depth = bcal_moth_Rio_Z[1] - bcal_moth_Rio_Z[0];
1670  return true;
1671  }
1672  else{
1673  _DBG_<<"Wrong vector size for BCAL mother RioZ!!!"<<endl;
1674  bcal_depth = 0.0;
1675  return false;
1676  }
1677 
1678 }
1679 
1680 //---------------------------------
1681 // GetBCALPhiShift
1682 //---------------------------------
1683 bool DGeometry::GetBCALPhiShift(float &bcal_phi_shift) const
1684 {
1685  vector<float> Phi0;
1686  bool good = Get("//section[@name='BarrelEMcal']/composition/mposPhi/@Phi0", Phi0);
1687  if(!good) return false;
1688  if(Phi0.size() == 1){
1689  bcal_phi_shift = Phi0[0];
1690  return true;
1691  }else{
1692  bcal_phi_shift = 0.0;
1693  return false;
1694  }
1695 }
1696 
1697 //---------------------------------
1698 // GetCCALZ
1699 //---------------------------------
1700 bool DGeometry::GetCCALZ(double &z_ccal) const
1701 {
1702  vector<double> ComptonEMcalpos;
1703  bool good = Get("//section/composition/posXYZ[@volume='ComptonEMcal']/@X_Y_Z", ComptonEMcalpos);
1704 
1705  if(!good){
1706  _DBG_<<"Unable to retrieve ComptonEMcal position."<<endl;
1707  z_ccal=876.106; // from some version of HDDS that may be out of date. 2018-12-10 DL
1708  return false;
1709  }else{
1710  z_ccal = ComptonEMcalpos[2];
1711  return true;
1712  }
1713 }
1714 
1715 //---------------------------------
1716 // GetFCALZ
1717 //---------------------------------
1718 bool DGeometry::GetFCALZ(double &z_fcal) const
1719 {
1720  vector<double> ForwardEMcalpos;
1721  bool good = Get("//section/composition/posXYZ[@volume='ForwardEMcal']/@X_Y_Z", ForwardEMcalpos);
1722 
1723  if(!good){
1724  _DBG_<<"Unable to retrieve ForwardEMcal position."<<endl;
1725  z_fcal=0.0;
1726  return false;
1727  }else{
1728  z_fcal = ForwardEMcalpos[2];
1729  return true;
1730  }
1731 }
1732 //---------------------------------
1733 // GetDIRCZ
1734 //---------------------------------
1735 bool DGeometry::GetDIRCZ(double &z_dirc) const
1736 {
1737  vector<double> dirc_face;
1738  bool good = Get("//section/composition/posXYZ[@volume='DIRC']/@X_Y_Z",dirc_face);
1739 
1740  if(!good){
1741  _DBG_<<"Unable to retrieve DIRC position."<<endl;
1742  z_dirc=0.0;
1743  return false;
1744  }
1745  else{
1746  vector<double>dirc_plane;
1747  vector<double>dirc_shift;
1748  vector<double>bar_plane;
1749  Get("//composition[@name='DRCC']/posXYZ[@volume='DCML10']/@X_Y_Z/plane[@value='1']", dirc_plane);
1750  Get("//composition[@name='DIRC']/posXYZ[@volume='DRCC']/@X_Y_Z", dirc_shift);
1751  z_dirc=dirc_face[2]+dirc_plane[2]+dirc_shift[2] + 0.8625; // last shift is the average center of quartz bar (585.862)
1752 
1753  jout << "DIRC z position = " << z_dirc << " cm." << endl;
1754  return true;
1755  }
1756 }
1757 
1758 //---------------------------------
1759 // GetTOFZ
1760 //---------------------------------
1761 bool DGeometry::GetTOFZ(vector<double> &z_tof) const
1762 {
1763  vector<double> ForwardTOF;
1764  vector<double> forwardTOF[2];
1765  vector<double> FTOC;
1766 
1767  if(!Get("//section/composition/posXYZ[@volume='ForwardTOF']/@X_Y_Z", ForwardTOF)) return false;
1768  if(!Get("//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='0']", forwardTOF[0])) return false;
1769  if(!Get("//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='1']", forwardTOF[1])) return false;
1770  if(!Get("//box[@name='FTOC' and sensitive='true']/@X_Y_Z", FTOC)) return false;
1771 
1772  z_tof.push_back(ForwardTOF[2] + forwardTOF[0][2] - FTOC[2]/2.0);
1773  z_tof.push_back(ForwardTOF[2] + forwardTOF[1][2] - FTOC[2]/2.0);
1774 
1775  //cerr << "DGeometry::GetTOFZ() = " << z_tof[0] << " " << z_tof[1] << endl;
1776 
1777  return true;
1778 }
1779 
1780 //---------------------------------
1781 // GetTOFPaddleParameters
1782 //---------------------------------
1783 bool DGeometry::GetTOFPaddleParameters(map<string,double> &paddle_params) const
1784 {
1785  vector<double> xyz_bar;
1786 
1787  // load the number of bars in each area
1788  int num_bars1 = 0;
1789  if(!Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOC']/@ncopy",num_bars1)) return false;
1790  int num_narrow_bars1 = 0;
1791  if(!Get("//composition[@name='forwardTOF_bottom2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars1)) return false;
1792  int num_single_end_bars1 = 0;
1793  if(!Get("//composition[@name='forwardTOF_north']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars1)) return false;
1794  int num_narrow_bars2 = 0;
1795  if(!Get("//composition[@name='forwardTOF_top2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars2)) return false;
1796  int num_bars2 = 0;
1797  if(!Get("//composition[@name='forwardTOF_top1']/mposY[@volume='FTOC']/@ncopy",num_bars2)) return false;
1798  int num_single_end_bars2 = 0;
1799  if(!Get("//composition[@name='forwardTOF_south']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars2)) return false;
1800 
1801  int NLONGBARS = num_bars1 + num_bars2 + num_narrow_bars1 + num_narrow_bars2;
1802  int NSHORTBARS = num_single_end_bars1 + num_single_end_bars2;
1803  int FIRSTSHORTBAR = num_bars1 + num_narrow_bars1 + 1;
1804  int LASTSHORTBAR = FIRSTSHORTBAR + NSHORTBARS/4;
1805 
1806  // load bar sizes
1807  //Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOC']/@X_Y_Z",xyz_bar);
1808  if(!Get("//box[@name='FTOC' and sensitive='true']/@X_Y_Z", xyz_bar)) return false;
1809  double LONGBARLENGTH = xyz_bar[0];
1810  double BARWIDTH = xyz_bar[1];
1811  //Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOH']/@X_Y_Z",xyz_bar);
1812  if(!Get("//box[@name='FTOH' and sensitive='true']/@X_Y_Z", xyz_bar)) return false;
1813  double SHORTBARLENGTH = xyz_bar[0];
1814 
1815 
1816  // load up the structure containing the parameters for the calling function
1817  paddle_params["NLONGBARS"] = NLONGBARS;
1818  paddle_params["NSHORTBARS"] = NSHORTBARS;
1819  paddle_params["BARWIDTH"] = BARWIDTH;
1820 
1821  paddle_params["LONGBARLENGTH"] = LONGBARLENGTH;
1822  paddle_params["HALFLONGBARLENGTH"] = LONGBARLENGTH/2.;
1823  paddle_params["SHORTBARLENGTH"] = SHORTBARLENGTH;
1824  paddle_params["HALFSHORTBARLENGTH"] = SHORTBARLENGTH/2.;
1825 
1826  paddle_params["FIRSTSHORTBAR"] = FIRSTSHORTBAR;
1827  paddle_params["LASTSHORTBAR"] = LASTSHORTBAR;
1828 
1829  //cout << "In DGeometry::GetTOFPaddleParameters() ..." << endl;
1830  //for(auto el : paddle_params) {
1831  // std::cout << el.first << " " << el.second << endl;
1832  //}
1833 
1834  return true;
1835 }
1836 
1837 
1838 //---------------------------------
1839 // GetTOFPaddlePerpPositions
1840 //---------------------------------
1841 bool DGeometry::GetTOFPaddlePerpPositions(vector<double> &y_tof) const
1842 {
1843  // add in a dummy entry, since we are indexing by paddle number, which starts at 1
1844  // maybe change this some day?
1845  y_tof.push_back(0);
1846 
1847  // Next fill array of bar positions within a plane
1848  // y_tof[barnumber] gives y position in the center of the bar. [currently barnumber = 1 - 46]
1849  double y0,dy;
1850 
1851  // load the number of bars
1852  int num_bars=1; // start counting at 1
1853  int num_bars1 = 0;
1854  Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOC']/@ncopy",num_bars1);
1855  int num_narrow_bars1 = 0;
1856  Get("//composition[@name='forwardTOF_bottom2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars1);
1857  int num_single_end_bars1 = 0;
1858  Get("//composition[@name='forwardTOF_north']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars1);
1859  int num_narrow_bars2 = 0;
1860  Get("//composition[@name='forwardTOF_top2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars2);
1861  int num_bars2 = 0;
1862  Get("//composition[@name='forwardTOF_top1']/mposY[@volume='FTOC']/@ncopy",num_bars2);
1863  int num_single_end_bars2 = 0;
1864  Get("//composition[@name='forwardTOF_south']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars2);
1865 
1866  // First 19 long bars
1867  Get("//composition[@name='forwardTOF_bottom1']/mposY/@Y0",y0);
1868  Get("//composition[@name='forwardTOF_bottom1']/mposY/@dY",dy);
1869  vector<double>tof_bottom1;
1870  Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_bottom1']/@X_Y_Z",tof_bottom1);
1871  for (int k=num_bars;k<num_bars+num_bars1;k++){
1872  y_tof.push_back(y0+tof_bottom1[1]+dy*double(k-1));
1873  }
1874  num_bars+=num_bars1;
1875 
1876  // two narrow long bars
1877  Get("//composition[@name='forwardTOF_bottom2']/mposY/@Y0",y0);
1878  Get("//composition[@name='forwardTOF_bottom2']/mposY/@dY",dy);
1879  vector<double>tof_bottom2;
1880  Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_bottom2']/@X_Y_Z",tof_bottom2);
1881  for (int k=num_bars;k<num_bars+num_narrow_bars1;k++){
1882  y_tof.push_back(y0+tof_bottom2[1]+dy*double(k-20));
1883  }
1884  num_bars+=num_narrow_bars1;
1885 
1886  // two short wide bars
1887  Get("//composition[@name='forwardTOF_north']/mposY/@Y0",y0);
1888  Get("//composition[@name='forwardTOF_north']/mposY/@dY",dy);
1889  vector<double>tof_north;
1890  Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_north']/@X_Y_Z",tof_north);
1891  for (int k=num_bars;k<num_bars+num_single_end_bars1;k++){
1892  y_tof.push_back(y0+tof_north[1]+dy*double(k-22));
1893  }
1894  num_bars+=num_single_end_bars1;
1895 
1896  // two narrow long bars
1897  Get("//composition[@name='forwardTOF_top2']/mposY/@Y0",y0);
1898  Get("//composition[@name='forwardTOF_top2']/mposY/@dY",dy);
1899  vector<double>tof_top2;
1900  Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_top2']/@X_Y_Z",tof_top2);
1901  for (int k=num_bars;k<num_bars+num_narrow_bars2;k++){
1902  y_tof.push_back(y0+tof_top2[1]+dy*double(k-24));
1903  }
1904  num_bars+=num_narrow_bars2;
1905 
1906  // Last 19 long bars
1907  Get("//composition[@name='forwardTOF_top1']/mposY/@Y0",y0);
1908  Get("//composition[@name='forwardTOF_top1']/mposY/@dY",dy);
1909  vector<double>tof_top1;
1910  Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_top1']/@X_Y_Z",tof_top1);
1911  for (int k=num_bars;k<num_bars+num_bars2;k++){
1912  y_tof.push_back(y0+tof_top1[1]+dy*double(k-26));
1913  }
1914  num_bars+=num_bars2;
1915 
1916  /*
1917  // two more short wide bars - IGNORE FOR NOW, ASSUME SAME Y AS OTHER SINGLE ENDED
1918  Get("//composition[@name='forwardTOF_south']/mposY/@Y0",y0);
1919  Get("//composition[@name='forwardTOF_south']/mposY/@dY",dy);
1920  vector<double>tof_south;
1921  Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_south']/@X_Y_Z",tof_south);
1922  for (unsigned int k=45;k<47;k++){
1923  y_tof.push_back(y0+tof_south[1]+dy*double(k-45));
1924  }
1925  */
1926 
1927  return true;
1928 }
1929 
1930 //---------------------------------
1931 // GetTargetZ
1932 //---------------------------------
1933 bool DGeometry::GetTargetZ(double &z_target) const
1934 {
1935  // Default to nominal center of GlueX target
1936  z_target=65.;
1937 
1938  // Check GlueX target is defined
1939  bool gluex_target_exists = true;
1940  vector<double> xyz_vessel;
1941  vector<double> xyz_target;
1942  vector<double> xyz_detector;
1943  if(gluex_target_exists) gluex_target_exists = Get("//composition[@name='targetVessel']/posXYZ[@volume='targetTube']/@X_Y_Z", xyz_vessel);
1944  if(gluex_target_exists) gluex_target_exists = Get("//composition[@name='Target']/posXYZ[@volume='targetVessel']/@X_Y_Z", xyz_target);
1945  if(gluex_target_exists) gluex_target_exists = Get("//posXYZ[@volume='Target']/@X_Y_Z", xyz_detector);
1946  if(gluex_target_exists) {
1947  z_target = xyz_vessel[2] + xyz_target[2] + xyz_detector[2];
1948  return true;
1949  }
1950 
1951  // Check if CPP target is defined
1952  bool cpp_target_exists = true;
1953  vector<double> xyz_TGT0;
1954  vector<double> xyz_TARG;
1955  vector<double> xyz_TargetCPP;
1956  if(cpp_target_exists) cpp_target_exists = Get("//composition/posXYZ[@volume='TGT0']/@X_Y_Z", xyz_TGT0);
1957  if(cpp_target_exists) cpp_target_exists = Get("//composition/posXYZ[@volume='TARG']/@X_Y_Z", xyz_TARG);
1958  if(cpp_target_exists) cpp_target_exists = Get("//composition/posXYZ[@volume='TargetCPP']/@X_Y_Z", xyz_TargetCPP);
1959  if(cpp_target_exists) {
1960  z_target = xyz_TGT0[2] + xyz_TARG[2] + xyz_TargetCPP[2];
1961  return true;
1962  }
1963 
1964  jout << " WARNING: Unable to get target location from XML for any of GlueX, or CPP targets. Using default of " << z_target << " cm" << endl;
1965 
1966  return false;
1967 }
1968 
1969 //---------------------------------
1970 // GetTargetLength
1971 //---------------------------------
1972 bool DGeometry::GetTargetLength(double &target_length) const
1973 {
1974  vector<double> zLength;
1975  bool good = Get("//section[@name='Target']/pcon[@name='LIH2']/real[@name='length']/[@value]", zLength);
1976 
1977  target_length = good ? zLength[0]:0.0;
1978 
1979  return good;
1980 }
1981 
1982 // Get vectors of positions and norm vectors for start counter from XML
1983 bool DGeometry::GetStartCounterGeom(vector<vector<DVector3> >&pos,
1984  vector<vector<DVector3> >&norm
1985  ) const{
1986 
1987  JCalibration *jcalib = dapp->GetJCalibration(runnumber);
1988 
1989  // Check if Start Counter geometry is present
1990  vector<double> sc_origin;
1991  bool got_sc = Get("//posXYZ[@volume='StartCntr']/@X_Y_Z", sc_origin);
1992  if(got_sc){
1993  // Offset from beam center
1994  vector<double>sc_origin_delta;
1995  Get("//posXYZ[@volume='startCntr']/@X_Y_Z", sc_origin_delta);
1996  double dx=sc_origin_delta[0];
1997  double dy=sc_origin_delta[1];
1998 
1999  // z-position at upstream face of scintillators.
2000  double z0=sc_origin[2];
2001 
2002  // Get rotation angles
2003  vector<double>sc_rot_angles;
2004  Get("//posXYZ[@volume='startCntr']/@rot", sc_rot_angles);
2005  double ThetaX=sc_rot_angles[0]*M_PI/180.;
2006  double ThetaY=sc_rot_angles[1]*M_PI/180.;
2007  Get("//posXYZ[@volume='StartCntr']/@rot", sc_rot_angles);
2008  //double ThetaX=sc_rot_angles[0]*M_PI/180.;
2009  //double ThetaY=sc_rot_angles[1]*M_PI/180.;
2010  double ThetaZ=sc_rot_angles[2]*M_PI/180.;
2011 
2012  // Get overall alignment shifts from CCDB
2013  map<string,double> sc_global_offsets;
2014  if (jcalib->Get("START_COUNTER/global_alignment_parms",sc_global_offsets)==false) {
2015  // translations
2016  dx += sc_global_offsets["SC_ALIGN_X"];
2017  dy += sc_global_offsets["SC_ALIGN_Y"];
2018  z0 += sc_global_offsets["SC_ALIGN_Z"];
2019 
2020  // rotations
2021  ThetaX += sc_global_offsets["SC_ALIGN_ROTX"]*M_PI/180.;
2022  ThetaY += sc_global_offsets["SC_ALIGN_ROTY"]*M_PI/180.;
2023  ThetaZ += sc_global_offsets["SC_ALIGN_ROTZ"]*M_PI/180.;
2024  }
2025 
2026  double num_paddles;
2027  Get("//mposPhi[@volume='STRC']/@ncopy",num_paddles);
2028  double dSCdphi = M_TWO_PI/num_paddles;
2029 
2030  vector<vector<double> > sc_rioz;
2031  GetMultiple("//pgon[@name='STRC']/polyplane/@Rio_Z", sc_rioz);
2032 
2033  // Get individual paddle alignment parameters
2034  vector< map<string,double> > sc_paddle_offsets;
2035  bool loaded_paddle_offsets = false;
2036  if (jcalib->Get("START_COUNTER/paddle_alignment_parms",sc_paddle_offsets)==false)
2037  loaded_paddle_offsets = true;
2038 
2039  // Create vectors of positions and normal vectors for each paddle
2040  for (unsigned int i=0;i<30;i++){
2041  double phi=ThetaZ+dSCdphi*(double(i)+0.5);
2042  double sinphi=sin(phi);
2043  double cosphi=cos(phi);
2044  double r=0.5*(sc_rioz[0][0]+sc_rioz[0][1]);
2045  DVector3 oldray;
2046  // Rotate by phi and take into account the tilt
2047  double x=r*cosphi;//+dx;
2048  double y=r*sinphi;//+dy;
2049  DVector3 ray(x,y,sc_rioz[0][2]);
2050  ray.RotateX(ThetaX);
2051  ray.RotateY(ThetaY);
2052 
2053  // Create stl-vectors to store positions and norm vectors
2054  vector<DVector3>posvec;
2055  vector<DVector3>dirvec;
2056  // Loop over radial/z positions describing start counter geometry from xml
2057  for(unsigned int k = 1; k < sc_rioz.size(); ++k){
2058  oldray=ray;
2059  r=0.5*(sc_rioz[k][0]+sc_rioz[k][1]);
2060  // Point in plane of scintillator
2061  x=r*cosphi;//+dx;
2062  y=r*sinphi;//+dy;
2063  ray.SetXYZ(x,y,sc_rioz[k][2]);
2064  ray.RotateX(ThetaX);
2065  ray.RotateY(ThetaY);
2066  // Apply alignment parameters
2067  if(loaded_paddle_offsets) {
2068  // allow for a maximum extent of the paddle in z
2069  double max_z = sc_paddle_offsets[i]["SC_MAX_Z"];
2070  if(ray.Z() > max_z) {
2071  ray.SetZ(max_z);
2072  }
2073  // allow for a modification of the bend angle of the paddle (18.5 deg from horizontal)
2074  // this should just be a perturbation around this angle, so assume a linear interpolation
2075  double delta_theta = sc_paddle_offsets[i]["SC_CURVE_THETA"]; // in degrees, r of curvature = 120 cm
2076  ray.SetX(ray.X()+delta_theta*1.65); // 1 degree ~ 1.65 cm in x
2077  ray.SetY(ray.Y()-delta_theta*0.55); // 1 degree ~ 0.55 cm in y
2078  }
2079  // Second point in the plane of the scintillator
2080  DVector3 ray2(r,10.,sc_rioz[k][2]);
2081  ray2.RotateZ(phi+0.5*dSCdphi*(1.+1./15.*((i>14)?29-i:i)));
2082  ray2.RotateX(ThetaX);
2083  ray2.RotateY(ThetaY);
2084  // Compute normal vector to plane
2085  DVector3 dir=(ray-oldray).Cross(ray2-oldray);
2086  dir.SetMag(1.);
2087  dirvec.push_back(dir);
2088  posvec.push_back(DVector3(oldray.X()+dx,oldray.Y()+dy,oldray.Z()+z0));
2089  }
2090  posvec.push_back(DVector3(ray.X(),ray.Y(),ray.Z()+z0)); //SAVE THE ENDPOINT OF THE LAST PLANE
2091  pos.push_back(posvec);
2092  norm.push_back(dirvec);
2093 
2094  posvec.clear();
2095  dirvec.clear();
2096  }
2097 
2098  }
2099  return got_sc;
2100 }
2101 
DApplication * dapp
bool GetFDCZ(vector< double > &z_wires) const
z-locations for each of the FDC wire planes in cm
Definition: DGeometry.cc:1221
double stereo_raw
Definition: DCDCWire.h:88
double rmax
#define X0
bool GetCCALZ(double &z_ccal) const
Definition: DGeometry.cc:1700
string toString(void)
Definition: DMaterial.cc:32
void FindNodes(string xpath, vector< xpathparsed_t > &matched_xpaths) const
Definition: DGeometry.cc:176
DVector3 angles
radians
Definition: DFDCWire.h:20
double phiZ
Definition: DCDCWire.h:85
float angle
radians
Definition: DFDCWire.h:19
TVector3 DVector3
Definition: DVector3.h:14
bool GetBCALRmin(float &bcal_rmin) const
minimum distance of BCAL module from beam line
Definition: DGeometry.cc:1544
Double_t x[NCHANNELS]
Definition: st_tw_resols.C:39
vector< DMaterialMap * > GetMaterialMapVector(void) const
Definition: DGeometry.cc:81
bool GetCompositeMaterial(const string &name, double &density, double &radlen) const
Definition: DGeometry.cc:527
double GetNr(void) const
Definition: DMaterialMap.h:59
sprintf(text,"Post KinFit Cut")
#define c
#define y
static vector< vector< DFDCWire * > > fdcwires
DMagneticFieldMap * GetBfield(unsigned int run_number=1)
bool GetFDCRmin(vector< double > &rmin_packages) const
beam hole size for each FDC package in cm
Definition: DGeometry.cc:1361
#define WIRES_PER_PLANE
Definition: DFDCGeometry.h:23
jerror_t FindMat(DVector3 &pos, double &rhoZ_overA, double &rhoZ_overA_logI, double &RadLen) const
Definition: DGeometry.cc:327
double zmax
int offsets(TString dir)
bool GetCDCOption(string &cdc_option) const
get the centralDC_option-X string
Definition: DGeometry.cc:1402
DLorentzDeflections * GetLorentzDeflections(unsigned int run_number=1)
double z0
Definition: DCDCWire.h:86
double rmin
bool GetFCALZ(double &z_fcal) const
z-location of front face of CCAL in cm
Definition: DGeometry.cc:1718
float angle
radians
Definition: DFDCCathode.h:15
bool GetFDCRmax(double &rmax_active_fdc) const
outer radius of FDC active area in cm
Definition: DGeometry.cc:1381
bool GetBCALDepth(float &bcal_depth) const
depth (or height) of BCAL module in cm
Definition: DGeometry.cc:1659
bool GetCDCWires(vector< vector< DCDCWire * > > &cdcwires) const
Definition: DGeometry.cc:773
int layer
1-24
Definition: DFDCWire.h:16
bool GetCDCAxialLength(double &cdc_axial_length) const
length of CDC axial wires in cm
Definition: DGeometry.cc:1425
DMagneticFieldMap * GetBfield(void) const
Definition: DGeometry.cc:61
void GetMaterials(void) const
Definition: DGeometry.cc:420
double r0
Definition: DCDCWire.h:89
bool GetTOFPaddlePerpPositions(vector< double > &y_tof) const
Definition: DGeometry.cc:1841
#define STRIP_SPACING
Definition: DFDCGeometry.h:27
bool GetBCALLength(float &bcal_length) const
length of BCAL module in cm
Definition: DGeometry.cc:1635
int straw
Definition: DCDCWire.h:81
bool GetBCALCenterZ(float &bcal_center_z) const
z-location of center of BCAL module in cm
Definition: DGeometry.cc:1611
double udir_mag
Definition: DCDCWire.h:87
jerror_t FindMatALT1(DVector3 &pos, DVector3 &mom, double &KrhoZ_overA, double &rhoZ_overA, double &LnI, double &X0, double *s_to_boundary=NULL) const
Definition: DGeometry.cc:221
bool GetBCALNmodules(unsigned int &bcal_nmodules) const
Number of BCAL modules.
Definition: DGeometry.cc:1589
virtual ~DGeometry()
Definition: DGeometry.cc:45
#define _DBG_
Definition: HDEVIO.h:12
bool GetCDCStereo(vector< double > &cdc_stereo) const
stereo angle for each CDC layer in degrees
Definition: DGeometry.cc:1441
double phiY
Definition: DCDCWire.h:85
#define FDC_ACTIVE_RADIUS
Definition: DFDCGeometry.h:14
#define U_OF_WIRE_ZERO
Definition: DFDCGeometry.h:25
void ReadMaterialMaps(void) const
Definition: DGeometry.cc:91
double phiX
Definition: DCDCWire.h:85
bool GetTOFZ(vector< double > &z_tof) const
z-location of front face of each of TOF in cm
Definition: DGeometry.cc:1761
double GetNz(void) const
Definition: DMaterialMap.h:60
double sqrt(double)
double sin(double)
bool GetBCALPhiShift(float &bcal_phi_shift) const
phi angle in degrees that first BCAL module is shifted from being centered at ph=0.0
Definition: DGeometry.cc:1683
bool GetTOFPaddleParameters(map< string, double > &paddle_params) const
Definition: DGeometry.cc:1783
#define FDC_NUM_LAYERS
Definition: DFDCGeometry.h:12
vector< node_t > xpathparsed_t
Definition: DGeometry.h:82
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
Extract the axial wire data from the XML.
Definition: DGeometry.cc:688
The DMaterial class holds information on a single material type. The main purpose is to hold informat...
Definition: DMaterial.h:25
bool IsInMap(const DVector3 &pos) const
bool GetFDCCathodes(vector< vector< DFDCCathode * > > &fdccathodes) const
Definition: DGeometry.cc:978
bool GetCDCNwires(vector< int > &cdc_nwires) const
Number of wires for each CDC layer.
Definition: DGeometry.cc:1459
float u
Definition: DFDCWire.h:18
#define M_TWO_PI
Definition: DGeometry.cc:20
int wire
1-N
Definition: DFDCWire.h:17
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
Find the materials traversed by a particle swimming from a specific position with a specific momentum...
Definition: DGeometry.cc:583
int ring
Definition: DCDCWire.h:80
#define STRIPS_PER_PLANE
Definition: DFDCGeometry.h:26
int layer
1-48
Definition: DFDCCathode.h:12
int strip
1-N
Definition: DFDCCathode.h:13
bool GetFDCWires(vector< vector< DFDCWire * > > &fdcwires) const
Definition: DGeometry.cc:1078
DLorentzDeflections * GetLorentzDeflections(void)
Definition: DGeometry.cc:73
bool GetCDCRmid(vector< double > &cdc_rmid) const
Distance of the center of CDC wire from beamline for each layer in cm.
Definition: DGeometry.cc:1450
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) const
Definition: DGeometry.cc:254
bool GetFDCStereo(vector< double > &stereo_angles) const
stereo angles of each of the FDC wire layers
Definition: DGeometry.cc:1299
bool GetCDCEndplate(double &z, double &dz, double &rmin, double &rmax) const
Definition: DGeometry.cc:1497
bool GetDIRCZ(double &z_dirc) const
z-location of DIRC in cm
Definition: DGeometry.cc:1735
bool GetBCALfADCRadii(vector< float > &fADC_radii) const
fADC radii including the outer radius of the last layer
Definition: DGeometry.cc:1567
#define CATHODE_ROT_ANGLE
Definition: DFDCGeometry.h:18
double y0
Definition: DCDCWire.h:86
TDirectory * dir
Definition: bcal_hist_eff.C:31
bool GetTargetLength(double &target_length) const
z-location of center of target
Definition: DGeometry.cc:1972
float phi
Definition: DCDCWire.h:78
double zmin
TH2F * temp
bool GetTargetZ(double &z_target) const
z-location of center of target
Definition: DGeometry.cc:1933
union @6::@8 u
#define WIRE_SPACING
Definition: DFDCGeometry.h:24
float stereo
Definition: DCDCWire.h:82
bool GetCDCCenterZ(double &cdc_center_z) const
z-location of center of CDC wires in cm
Definition: DGeometry.cc:1416
bool GetStartCounterGeom(vector< vector< DVector3 > > &pos, vector< vector< DVector3 > > &norm) const
Definition: DGeometry.cc:1983
const DMaterialMap * FindDMaterialMap(DVector3 &pos) const
Definition: DGeometry.cc:363
double x0
Definition: DCDCWire.h:86
const DMaterial * GetDMaterial(string name) const
Definition: DGeometry.cc:396
double phiStraw
Definition: DCDCWire.h:90