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