File: | /home/sdobbs/work/clang/halld_recon/src/libraries/HDGEOMETRY/DGeometry.cc |
Warning: | line 2116, column 40 The left operand of '*' is a garbage 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)found_all = GetCompositeMaterial(name, density, radlen); | |||
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)continue; | |||
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_ |