| File: | /volatile/halld/gluex/nightly/2024-03-18/Linux_CentOS7.7-x86_64-gcc4.8.5/halld_recon/src/libraries/HDGEOMETRY/DGeometry.cc |
| Location: | line 2356, column 4 |
| Description: | Value stored to 'num_bars' is never read |
| 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 | |
| 1788 | //--------------------------------- |
| 1789 | // GetCTOFZ |
| 1790 | //--------------------------------- |
| 1791 | bool DGeometry::GetCTOFZ(double &z) const { |
| 1792 | z=1000; // cm; initialize to a large value |
| 1793 | vector<double> CppScintPos; |
| 1794 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 1795 | bool good = Get("//section/composition/posXYZ[@volume='CppScint']/@X_Y_Z", CppScintPos); |
| 1796 | jgeom->SetVerbose(1); // reenable error messages |
| 1797 | if (!good){ |
| 1798 | //_DBG_<<"Unable to retrieve CPP scintillator position."<<endl; |
| 1799 | return false; |
| 1800 | } |
| 1801 | z=CppScintPos[2]; |
| 1802 | return true; |
| 1803 | } |
| 1804 | |
| 1805 | //--------------------------------- |
| 1806 | // GetCTOFPositions |
| 1807 | //--------------------------------- |
| 1808 | bool DGeometry::GetCTOFPositions(vector<DVector3>&posvec) const{ |
| 1809 | vector<double>origin; |
| 1810 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 1811 | bool good = Get("//section/composition/posXYZ[@volume='CppScint']/@X_Y_Z",origin); |
| 1812 | jgeom->SetVerbose(1); // reenable error messages |
| 1813 | if (!good) return false; |
| 1814 | DVector3 pos(origin[0],origin[1],origin[2]); |
| 1815 | for (unsigned int paddle=1;paddle<5;paddle++){ |
| 1816 | vector<double>local_pos; |
| 1817 | Get(Form("//posXYZ[@volume='CPPPaddle']/@X_Y_Z/column[@value='%d']",paddle),local_pos); |
| 1818 | DVector3 dpos(local_pos[0],local_pos[1],local_pos[2]); |
| 1819 | posvec.push_back(pos+dpos); |
| 1820 | } |
| 1821 | |
| 1822 | return true; |
| 1823 | } |
| 1824 | |
| 1825 | //--------------------------------- |
| 1826 | // GetFMWPCZ |
| 1827 | //--------------------------------- |
| 1828 | bool DGeometry::GetFMWPCZ_vec(vector<double>&zvec_fmwpc) const |
| 1829 | { |
| 1830 | vector<double> ForwardMWPCpos; |
| 1831 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 1832 | bool good = Get("//section/composition/posXYZ[@volume='ForwardMWPC']/@X_Y_Z", ForwardMWPCpos); |
| 1833 | jgeom->SetVerbose(1); // reenable error messages |
| 1834 | if (!good){ |
| 1835 | //_DBG_<<"Unable to retrieve ForwardMWPC position."<<endl; |
| 1836 | return false; |
| 1837 | } |
| 1838 | |
| 1839 | vector<double>CPPChamberPos; |
| 1840 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='1']", CPPChamberPos); |
| 1841 | zvec_fmwpc.push_back(ForwardMWPCpos[2]+CPPChamberPos[2]); |
| 1842 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='2']", CPPChamberPos); |
| 1843 | zvec_fmwpc.push_back(ForwardMWPCpos[2]+CPPChamberPos[2]); |
| 1844 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='3']", CPPChamberPos); |
| 1845 | zvec_fmwpc.push_back(ForwardMWPCpos[2]+CPPChamberPos[2]); |
| 1846 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='4']", CPPChamberPos); |
| 1847 | zvec_fmwpc.push_back(ForwardMWPCpos[2]+CPPChamberPos[2]); |
| 1848 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='5']", CPPChamberPos); |
| 1849 | zvec_fmwpc.push_back(ForwardMWPCpos[2]+CPPChamberPos[2]); |
| 1850 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='6']", CPPChamberPos); |
| 1851 | zvec_fmwpc.push_back(ForwardMWPCpos[2]+CPPChamberPos[2]); |
| 1852 | |
| 1853 | return true; |
| 1854 | } |
| 1855 | |
| 1856 | //--------------------------------- |
| 1857 | // GetFMWPCXY |
| 1858 | //--------------------------------- |
| 1859 | bool DGeometry::GetFMWPCXY_vec(vector<double>&xvec_fmwpc, vector<double>&yvec_fmwpc) const |
| 1860 | { |
| 1861 | vector<double> ForwardMWPCpos; |
| 1862 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 1863 | bool good = Get("//section/composition/posXYZ[@volume='ForwardMWPC']/@X_Y_Z", ForwardMWPCpos); |
| 1864 | jgeom->SetVerbose(1); // reenable error messages |
| 1865 | if (!good){ |
| 1866 | //_DBG_<<"Unable to retrieve ForwardMWPC position."<<endl; |
| 1867 | return false; |
| 1868 | } |
| 1869 | |
| 1870 | // Get offsets tweaking nominal geometry from calibration database |
| 1871 | // JCalibration * jcalib = dapp->GetJCalibration(runnumber); |
| 1872 | |
| 1873 | vector<double>CPPChamberPos; |
| 1874 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='1']", CPPChamberPos); |
| 1875 | xvec_fmwpc.push_back(ForwardMWPCpos[0]+CPPChamberPos[0]); |
| 1876 | yvec_fmwpc.push_back(ForwardMWPCpos[1]+CPPChamberPos[1]); |
| 1877 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='2']", CPPChamberPos); |
| 1878 | xvec_fmwpc.push_back(ForwardMWPCpos[0]+CPPChamberPos[0]); |
| 1879 | yvec_fmwpc.push_back(ForwardMWPCpos[1]+CPPChamberPos[1]); |
| 1880 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='3']", CPPChamberPos); |
| 1881 | xvec_fmwpc.push_back(ForwardMWPCpos[0]+CPPChamberPos[0]); |
| 1882 | yvec_fmwpc.push_back(ForwardMWPCpos[1]+CPPChamberPos[1]); |
| 1883 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='4']", CPPChamberPos); |
| 1884 | xvec_fmwpc.push_back(ForwardMWPCpos[0]+CPPChamberPos[0]); |
| 1885 | yvec_fmwpc.push_back(ForwardMWPCpos[1]+CPPChamberPos[1]); |
| 1886 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='5']", CPPChamberPos); |
| 1887 | xvec_fmwpc.push_back(ForwardMWPCpos[0]+CPPChamberPos[0]); |
| 1888 | yvec_fmwpc.push_back(ForwardMWPCpos[1]+CPPChamberPos[1]); |
| 1889 | Get("//posXYZ[@volume='CPPChamber']/@X_Y_Z/layer[@value='6']", CPPChamberPos); |
| 1890 | xvec_fmwpc.push_back(ForwardMWPCpos[0]+CPPChamberPos[0]); |
| 1891 | yvec_fmwpc.push_back(ForwardMWPCpos[1]+CPPChamberPos[1]); |
| 1892 | |
| 1893 | // Currently, not all chambers have a 'rot' field in hdds |
| 1894 | // vector<double>CPPChamberRot; |
| 1895 | // Get("//posXYZ[@volume='CPPChamber']/@rot/layer[@value='1']", CPPChamberRot); |
| 1896 | // rot_fmwpc.push_back(CPPChamberRot[2]); |
| 1897 | // Get("//posXYZ[@volume='CPPChamber']/@rot/layer[@value='2']", CPPChamberRot); |
| 1898 | // rot_fmwpc.push_back(CPPChamberRot[2]); |
| 1899 | // Get("//posXYZ[@volume='CPPChamber']/@rot/layer[@value='3']", CPPChamberRot); |
| 1900 | // rot_fmwpc.push_back(CPPChamberRot[2]); |
| 1901 | // Get("//posXYZ[@volume='CPPChamber']/@rot/layer[@value='4']", CPPChamberRot); |
| 1902 | // rot_fmwpc.push_back(CPPChamberRot[2]); |
| 1903 | // Get("//posXYZ[@volume='CPPChamber']/@rot/layer[@value='5']", CPPChamberRot); |
| 1904 | // rot_fmwpc.push_back(CPPChamberRot[2]); |
| 1905 | // Get("//posXYZ[@volume='CPPChamber']/@rot/layer[@value='6']", CPPChamberRot); |
| 1906 | // rot_fmwpc.push_back(CPPChamberRot[2]); |
| 1907 | |
| 1908 | return true; |
| 1909 | } |
| 1910 | |
| 1911 | //--------------------------------- |
| 1912 | // GetFMWPCSize -- use the dimensions of the frame |
| 1913 | //--------------------------------- |
| 1914 | bool DGeometry::GetFMWPCSize(double &xy_fmwpc) const |
| 1915 | { |
| 1916 | vector<double> ForwardMWPCdimensions; |
| 1917 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 1918 | bool good = Get("//section[@name='ForwardMWPC']/box[@name='CPPF']/@X_Y_Z", ForwardMWPCdimensions); |
| 1919 | jgeom->SetVerbose(1); // reenable error messages |
| 1920 | if (!good){ |
| 1921 | xy_fmwpc=0.0; |
| 1922 | return false; |
| 1923 | } |
| 1924 | xy_fmwpc=0.5*ForwardMWPCdimensions[0]; |
| 1925 | |
| 1926 | return true; |
| 1927 | } |
| 1928 | |
| 1929 | //--------------------------------- |
| 1930 | // GetFMWPCWireSpacing -- space between wires in cm |
| 1931 | //--------------------------------- |
| 1932 | bool DGeometry::GetFMWPCWireSpacing(double &fmwpc_wire_spacing) const |
| 1933 | { |
| 1934 | fmwpc_wire_spacing = 1.016; |
| 1935 | |
| 1936 | return true; |
| 1937 | } |
| 1938 | |
| 1939 | //--------------------------------- |
| 1940 | // GetFMWPCWireSpacing -- space between wires in cm |
| 1941 | //--------------------------------- |
| 1942 | bool DGeometry::GetFMWPCWireOrientation(vector<fmwpc_wire_orientation_t> &fmwpc_wire_orientation) const |
| 1943 | { |
| 1944 | fmwpc_wire_orientation.clear(); |
| 1945 | fmwpc_wire_orientation.push_back( kFMWPC_WIRE_ORIENTATION_VERTICAL ); |
| 1946 | fmwpc_wire_orientation.push_back( kFMWPC_WIRE_ORIENTATION_HORIZONTAL ); |
| 1947 | fmwpc_wire_orientation.push_back( kFMWPC_WIRE_ORIENTATION_VERTICAL ); |
| 1948 | fmwpc_wire_orientation.push_back( kFMWPC_WIRE_ORIENTATION_HORIZONTAL ); |
| 1949 | fmwpc_wire_orientation.push_back( kFMWPC_WIRE_ORIENTATION_VERTICAL ); |
| 1950 | fmwpc_wire_orientation.push_back( kFMWPC_WIRE_ORIENTATION_HORIZONTAL ); |
| 1951 | |
| 1952 | return true; |
| 1953 | } |
| 1954 | |
| 1955 | //--------------------------------- |
| 1956 | // GetFCALZ |
| 1957 | //--------------------------------- |
| 1958 | bool DGeometry::GetFCALZ(double &z_fcal) const |
| 1959 | { |
| 1960 | vector<double> ForwardEMcalpos; |
| 1961 | bool good = Get("//section/composition/posXYZ[@volume='ForwardEMcal']/@X_Y_Z", ForwardEMcalpos); |
| 1962 | |
| 1963 | if(!good){ |
| 1964 | _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":" <<1964<<" "<<"Unable to retrieve ForwardEMcal position."<<endl; |
| 1965 | z_fcal=0.0; |
| 1966 | return false; |
| 1967 | }else{ |
| 1968 | z_fcal = ForwardEMcalpos[2]; |
| 1969 | return true; |
| 1970 | } |
| 1971 | } |
| 1972 | |
| 1973 | |
| 1974 | //--------------------------------- |
| 1975 | // GetFCALPosition |
| 1976 | //--------------------------------- |
| 1977 | bool DGeometry::GetFCALPosition(double &x,double &y,double &z) const |
| 1978 | { |
| 1979 | vector<double> ForwardEMcalpos; |
| 1980 | bool good = Get("//section/composition/posXYZ[@volume='ForwardEMcal']/@X_Y_Z", ForwardEMcalpos); |
| 1981 | |
| 1982 | if(!good){ |
| 1983 | _DBG_std::cerr<<"libraries/HDGEOMETRY/DGeometry.cc"<<":" <<1983<<" "<<"Unable to retrieve ForwardEMcal position."<<endl; |
| 1984 | x=0.,y=0.,z=0.; |
| 1985 | return false; |
| 1986 | }else{ |
| 1987 | x=ForwardEMcalpos[0],y=ForwardEMcalpos[1],z=ForwardEMcalpos[2]; |
| 1988 | //_DBG_ << "FCAL position: (x,y,z)=(" << x <<"," << y << "," << z << ")"<<endl; |
| 1989 | return true; |
| 1990 | } |
| 1991 | } |
| 1992 | |
| 1993 | //--------------------------------- |
| 1994 | // GetCCALPosition |
| 1995 | //--------------------------------- |
| 1996 | bool DGeometry::GetCCALPosition(double &x,double &y,double &z) const |
| 1997 | { |
| 1998 | vector<double> ComptonEMcalpos; |
| 1999 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 2000 | bool good = Get("//section/composition/posXYZ[@volume='ComptonEMcal']/@X_Y_Z", ComptonEMcalpos); |
| 2001 | jgeom->SetVerbose(1); // reenable error messages |
| 2002 | |
| 2003 | if(!good){ |
| 2004 | //_DBG_<<"Unable to retrieve ComptonEMcal position."<<endl; |
| 2005 | x=0.,y=0.,z=0.; |
| 2006 | return false; |
| 2007 | }else{ |
| 2008 | x=ComptonEMcalpos[0],y=ComptonEMcalpos[1],z=ComptonEMcalpos[2]; |
| 2009 | //_DBG_ << "CCAL position: (x,y,z)=(" << ComptonEMcalpos[0] <<"," |
| 2010 | // << ComptonEMcalpos[1]<<","<<ComptonEMcalpos[2]<< ")" << endl; |
| 2011 | return true; |
| 2012 | } |
| 2013 | } |
| 2014 | |
| 2015 | //--------------------------------- |
| 2016 | // GetFCALInsertRowSize |
| 2017 | //--------------------------------- |
| 2018 | bool DGeometry::GetFCALInsertRowSize(int &insert_row_size) const |
| 2019 | { |
| 2020 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 2021 | bool good = Get("//composition[@name='LeadTungstateFullRow']/mposX[@volume='LTBLwrapped']/@ncopy",insert_row_size); |
| 2022 | jgeom->SetVerbose(1); // reenable error messages |
| 2023 | |
| 2024 | if(!good){ |
| 2025 | // NEED TO RETHINK ERROR REPORTING FOR OPTIONAL DETECTOR ELEMENTS |
| 2026 | //_DBG_<<"Unable to retrieve ComptonEMcal position."<<endl; |
| 2027 | insert_row_size = 0; |
| 2028 | return false; |
| 2029 | }else{ |
| 2030 | return true; |
| 2031 | } |
| 2032 | } |
| 2033 | |
| 2034 | //--------------------------------- |
| 2035 | // GetFCALBlockSize |
| 2036 | //--------------------------------- |
| 2037 | bool DGeometry::GetFCALBlockSize(vector<double> &block) const |
| 2038 | { |
| 2039 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 2040 | bool good = Get("//box[@name='LGBL']/@X_Y_Z",block); |
| 2041 | jgeom->SetVerbose(1); // reenable error messages |
| 2042 | |
| 2043 | if(!good){ |
| 2044 | // NEED TO RETHINK ERROR REPORTING FOR OPTIONAL DETECTOR ELEMENTS |
| 2045 | //_DBG_<<"Unable to retrieve ComptonEMcal position."<<endl; |
| 2046 | return false; |
| 2047 | }else{ |
| 2048 | return true; |
| 2049 | } |
| 2050 | } |
| 2051 | |
| 2052 | //--------------------------------- |
| 2053 | // GetFCALInsertBlockSize |
| 2054 | //--------------------------------- |
| 2055 | bool DGeometry::GetFCALInsertBlockSize(vector<double> &block) const |
| 2056 | { |
| 2057 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 2058 | bool good = Get("//box[@name='LTB1']/@X_Y_Z",block); |
| 2059 | jgeom->SetVerbose(1); // reenable error messages |
| 2060 | |
| 2061 | if(!good){ |
| 2062 | // NEED TO RETHINK ERROR REPORTING FOR OPTIONAL DETECTOR ELEMENTS |
| 2063 | //_DBG_<<"Unable to retrieve ComptonEMcal position."<<endl; |
| 2064 | return false; |
| 2065 | }else{ |
| 2066 | return true; |
| 2067 | } |
| 2068 | } |
| 2069 | |
| 2070 | //--------------------------------- |
| 2071 | // GetDIRCZ |
| 2072 | //--------------------------------- |
| 2073 | bool DGeometry::GetDIRCZ(double &z_dirc) const |
| 2074 | { |
| 2075 | vector<double> dirc_face; |
| 2076 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 2077 | bool good = Get("//section/composition/posXYZ[@volume='DIRC']/@X_Y_Z",dirc_face); |
| 2078 | jgeom->SetVerbose(1); // reenable error messages |
| 2079 | |
| 2080 | if(!good){ |
| 2081 | //_DBG_<<"Unable to retrieve DIRC position."<<endl; |
| 2082 | z_dirc=0.0; |
| 2083 | return false; |
| 2084 | } |
| 2085 | else{ |
| 2086 | vector<double>dirc_plane; |
| 2087 | vector<double>dirc_shift; |
| 2088 | vector<double>bar_plane; |
| 2089 | Get("//composition[@name='DRCC']/posXYZ[@volume='DCML10']/@X_Y_Z/plane[@value='1']", dirc_plane); |
| 2090 | Get("//composition[@name='DIRC']/posXYZ[@volume='DRCC']/@X_Y_Z", dirc_shift); |
| 2091 | z_dirc=dirc_face[2]+dirc_plane[2]+dirc_shift[2] + 0.8625; // last shift is the average center of quartz bar (585.862) |
| 2092 | |
| 2093 | //jout << "DIRC z position = " << z_dirc << " cm." << endl; |
| 2094 | return true; |
| 2095 | } |
| 2096 | } |
| 2097 | |
| 2098 | //--------------------------------- |
| 2099 | // GetTRDZ |
| 2100 | //--------------------------------- |
| 2101 | bool DGeometry::GetTRDZ(vector<double> &z_trd) const |
| 2102 | { |
| 2103 | vector<double> trd_origin; |
| 2104 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 2105 | bool good = Get("//section/composition/posXYZ[@volume='TRDGEM']/@X_Y_Z",trd_origin); |
| 2106 | jgeom->SetVerbose(1); // reenable error messages |
| 2107 | |
| 2108 | if(!good){ |
| 2109 | //_DBG_<<"Unable to retrieve TRD position."<<endl; |
| 2110 | return false; |
| 2111 | } |
| 2112 | else{ |
| 2113 | vector<double>trd_G10; |
| 2114 | Get("//composition[@name='TRDGEM']/posXYZ[@volume='TGPK']/@X_Y_Z",trd_G10); |
| 2115 | |
| 2116 | jout << "TRD z positions = "; |
| 2117 | for(int i=0; i<5; i++) { |
| 2118 | vector<double>trd_plane; |
| 2119 | Get(Form("//composition[@name='TGPK']/posXYZ[@volume='TRDG']/@X_Y_Z/plane[@value='%d']",i),trd_plane); |
| 2120 | double z_trd_plane=trd_origin[2]+trd_G10[2]+trd_plane[2]; |
| 2121 | jout << z_trd_plane << ", "; |
| 2122 | z_trd.push_back(z_trd_plane); |
| 2123 | } |
| 2124 | jout << "cm" << endl; |
| 2125 | } |
| 2126 | |
| 2127 | return true; |
| 2128 | } |
| 2129 | //--------------------------------- |
| 2130 | // GetTOFZ |
| 2131 | //--------------------------------- |
| 2132 | bool DGeometry::GetTOFZ(double &CenterVPlane,double &CenterHPlane, |
| 2133 | double &CenterMPlane) const{ |
| 2134 | // Store the z position for both planes |
| 2135 | vector<double>tof_face; |
| 2136 | if (!Get("//section/composition/posXYZ[@volume='ForwardTOF']/@X_Y_Z",tof_face)){ |
| 2137 | return false; |
| 2138 | } |
| 2139 | vector<double>tof_plane0; |
| 2140 | Get("//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='0']", tof_plane0); |
| 2141 | vector<double>tof_plane1; |
| 2142 | Get("//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='1']", tof_plane1); |
| 2143 | CenterVPlane=tof_face[2]+tof_plane1[2]; |
| 2144 | CenterHPlane=tof_face[2]+tof_plane0[2]; |
| 2145 | // also save position midway between the two planes |
| 2146 | CenterMPlane=0.5*(CenterHPlane+CenterVPlane); |
| 2147 | |
| 2148 | return true; |
| 2149 | } |
| 2150 | |
| 2151 | //--------------------------------- |
| 2152 | // GetTOFZ |
| 2153 | //--------------------------------- |
| 2154 | bool DGeometry::GetTOFZ(vector<double> &z_tof) const |
| 2155 | { |
| 2156 | vector<double> ForwardTOF; |
| 2157 | vector<double> forwardTOF[2]; |
| 2158 | vector<double> FTOC; |
| 2159 | |
| 2160 | if(!Get("//section/composition/posXYZ[@volume='ForwardTOF']/@X_Y_Z", ForwardTOF)) return false; |
| 2161 | if(!Get("//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='0']", forwardTOF[0])) return false; |
| 2162 | if(!Get("//composition[@name='ForwardTOF']/posXYZ[@volume='forwardTOF']/@X_Y_Z/plane[@value='1']", forwardTOF[1])) return false; |
| 2163 | if(!Get("//box[@name='FTOC' and sensitive='true']/@X_Y_Z", FTOC)) return false; |
| 2164 | |
| 2165 | z_tof.push_back(ForwardTOF[2] + forwardTOF[0][2] - FTOC[2]/2.0); |
| 2166 | z_tof.push_back(ForwardTOF[2] + forwardTOF[1][2] - FTOC[2]/2.0); |
| 2167 | |
| 2168 | //cerr << "DGeometry::GetTOFZ() = " << z_tof[0] << " " << z_tof[1] << endl; |
| 2169 | |
| 2170 | return true; |
| 2171 | } |
| 2172 | |
| 2173 | //--------------------------------- |
| 2174 | // GetTOFPaddleParameters |
| 2175 | //--------------------------------- |
| 2176 | bool DGeometry::GetTOFPaddleParameters(map<string,double> &paddle_params) const |
| 2177 | { |
| 2178 | vector<double> xyz_bar; |
| 2179 | |
| 2180 | // load the number of bars in each area |
| 2181 | int num_bars1 = 0; |
| 2182 | if(!Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOC']/@ncopy",num_bars1)) return false; |
| 2183 | int num_narrow_bars1 = 0; |
| 2184 | if(!Get("//composition[@name='forwardTOF_bottom2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars1)) return false; |
| 2185 | int num_single_end_bars1 = 0; |
| 2186 | if(!Get("//composition[@name='forwardTOF_north']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars1)) return false; |
| 2187 | |
| 2188 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 2189 | int num_narrower_bars1 = 0; // optional - added during upgrade |
| 2190 | Get("//composition[@name='forwardTOF_bottom3']/mposY[@volume='FTOL']/@ncopy",num_narrower_bars1); |
| 2191 | int num_narrower_bars2 = 0; // optional - added during upgrade |
| 2192 | Get("//composition[@name='forwardTOF_top3']/mposY[@volume='FTOL']/@ncopy",num_narrower_bars2); |
| 2193 | jgeom->SetVerbose(1); // reenable error messages |
| 2194 | |
| 2195 | int num_narrow_bars2 = 0; |
| 2196 | if(!Get("//composition[@name='forwardTOF_top2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars2)) return false; |
| 2197 | int num_bars2 = 0; |
| 2198 | if(!Get("//composition[@name='forwardTOF_top1']/mposY[@volume='FTOC']/@ncopy",num_bars2)) return false; |
| 2199 | int num_single_end_bars2 = 0; |
| 2200 | if(!Get("//composition[@name='forwardTOF_south']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars2)) return false; |
| 2201 | |
| 2202 | int NLONGBARS = num_bars1 + num_bars2 + num_narrow_bars1 + num_narrow_bars2 |
| 2203 | + num_narrower_bars1 + num_narrower_bars2; |
| 2204 | int NSHORTBARS = num_single_end_bars1 + num_single_end_bars2; |
| 2205 | int FIRSTSHORTBAR = num_bars1 + num_narrow_bars1 + num_narrower_bars1 + 1; |
| 2206 | int LASTSHORTBAR = FIRSTSHORTBAR + NSHORTBARS/2 - 1; |
| 2207 | |
| 2208 | // load bar sizes |
| 2209 | //Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOC']/@X_Y_Z",xyz_bar); |
| 2210 | if(!Get("//box[@name='FTOC' and sensitive='true']/@X_Y_Z", xyz_bar)) return false; |
| 2211 | double LONGBARLENGTH = xyz_bar[0]; |
| 2212 | double BARWIDTH = xyz_bar[1]; |
| 2213 | //Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOH']/@X_Y_Z",xyz_bar); |
| 2214 | if(!Get("//box[@name='FTOH' and sensitive='true']/@X_Y_Z", xyz_bar)) return false; |
| 2215 | double SHORTBARLENGTH = xyz_bar[0]; |
| 2216 | |
| 2217 | |
| 2218 | // load up the structure containing the parameters for the calling function |
| 2219 | paddle_params["NLONGBARS"] = NLONGBARS; |
| 2220 | paddle_params["NSHORTBARS"] = NSHORTBARS; |
| 2221 | paddle_params["BARWIDTH"] = BARWIDTH; |
| 2222 | |
| 2223 | paddle_params["LONGBARLENGTH"] = LONGBARLENGTH; |
| 2224 | paddle_params["HALFLONGBARLENGTH"] = LONGBARLENGTH/2.; |
| 2225 | paddle_params["SHORTBARLENGTH"] = SHORTBARLENGTH; |
| 2226 | paddle_params["HALFSHORTBARLENGTH"] = SHORTBARLENGTH/2.; |
| 2227 | |
| 2228 | paddle_params["FIRSTSHORTBAR"] = FIRSTSHORTBAR; |
| 2229 | paddle_params["LASTSHORTBAR"] = LASTSHORTBAR; |
| 2230 | |
| 2231 | //cout << "In DGeometry::GetTOFPaddleParameters() ..." << endl; |
| 2232 | //for(auto el : paddle_params) { |
| 2233 | // std::cout << el.first << " " << el.second << endl; |
| 2234 | //} |
| 2235 | |
| 2236 | return true; |
| 2237 | } |
| 2238 | |
| 2239 | |
| 2240 | //--------------------------------- |
| 2241 | // GetTOFPaddlePerpPositions |
| 2242 | //--------------------------------- |
| 2243 | bool DGeometry::GetTOFPaddlePerpPositions(vector<double> &y_tof, vector<double> &y_widths) const |
| 2244 | { |
| 2245 | // add in a dummy entry, since we are indexing by paddle number, which starts at 1 |
| 2246 | // maybe change this some day? |
| 2247 | y_tof.push_back(0); |
| 2248 | y_widths.push_back(0); |
| 2249 | |
| 2250 | // Next fill array of bar positions within a plane |
| 2251 | // y_tof[barnumber] gives y position in the center of the bar. [currently barnumber = 1 - 46] |
| 2252 | double y0,dy; |
| 2253 | |
| 2254 | // load the number of bars |
| 2255 | int num_bars=1; // start counting at 1 |
| 2256 | int num_bars1 = 0; |
| 2257 | Get("//composition[@name='forwardTOF_bottom1']/mposY[@volume='FTOC']/@ncopy",num_bars1); |
| 2258 | int num_narrow_bars1 = 0; |
| 2259 | Get("//composition[@name='forwardTOF_bottom2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars1); |
| 2260 | int num_single_end_bars1 = 0; |
| 2261 | Get("//composition[@name='forwardTOF_north']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars1); |
| 2262 | |
| 2263 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 2264 | int num_narrower_bars1 = 0; |
| 2265 | Get("//composition[@name='forwardTOF_bottom3']/mposY[@volume='FTOL']/@ncopy",num_narrower_bars1); |
| 2266 | int num_narrower_bars2 = 0; |
| 2267 | Get("//composition[@name='forwardTOF_top3']/mposY[@volume='FTOL']/@ncopy",num_narrower_bars2); |
| 2268 | jgeom->SetVerbose(1); // reenable error messages |
| 2269 | |
| 2270 | int num_narrow_bars2 = 0; |
| 2271 | Get("//composition[@name='forwardTOF_top2']/mposY[@volume='FTOX']/@ncopy",num_narrow_bars2); |
| 2272 | int num_bars2 = 0; |
| 2273 | Get("//composition[@name='forwardTOF_top1']/mposY[@volume='FTOC']/@ncopy",num_bars2); |
| 2274 | int num_single_end_bars2 = 0; |
| 2275 | Get("//composition[@name='forwardTOF_south']/mposY[@volume='FTOH']/@ncopy",num_single_end_bars2); |
| 2276 | |
| 2277 | // First 19 long bars |
| 2278 | Get("//composition[@name='forwardTOF_bottom1']/mposY/@Y0",y0); |
| 2279 | Get("//composition[@name='forwardTOF_bottom1']/mposY/@dY",dy); |
| 2280 | vector<double>tof_bottom1; |
| 2281 | Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_bottom1']/@X_Y_Z",tof_bottom1); |
| 2282 | for (int k=num_bars;k<num_bars+num_bars1;k++){ |
| 2283 | y_tof.push_back(y0+tof_bottom1[1]+dy*double(k-num_bars)); |
| 2284 | y_widths.push_back(dy); |
| 2285 | } |
| 2286 | num_bars+=num_bars1; |
| 2287 | |
| 2288 | // two narrow long bars |
| 2289 | Get("//composition[@name='forwardTOF_bottom2']/mposY/@Y0",y0); |
| 2290 | Get("//composition[@name='forwardTOF_bottom2']/mposY/@dY",dy); |
| 2291 | vector<double>tof_bottom2; |
| 2292 | Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_bottom2']/@X_Y_Z",tof_bottom2); |
| 2293 | for (int k=num_bars;k<num_bars+num_narrow_bars1;k++){ |
| 2294 | y_tof.push_back(y0+tof_bottom2[1]+dy*double(k-num_bars)); |
| 2295 | y_widths.push_back(dy); |
| 2296 | } |
| 2297 | num_bars+=num_narrow_bars1; |
| 2298 | |
| 2299 | // two narrower long bars - added by upgrade |
| 2300 | if(num_narrower_bars1 > 0) { |
| 2301 | Get("//composition[@name='forwardTOF_bottom3']/mposY/@Y0",y0); |
| 2302 | Get("//composition[@name='forwardTOF_bottom3']/mposY/@dY",dy); |
| 2303 | vector<double>tof_bottom3; |
| 2304 | Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_bottom3']/@X_Y_Z",tof_bottom3); |
| 2305 | for (int k=num_bars;k<num_bars+num_narrower_bars1;k++){ |
| 2306 | y_tof.push_back(y0+tof_bottom3[1]+dy*double(k-num_bars)); |
| 2307 | y_widths.push_back(dy); |
| 2308 | } |
| 2309 | num_bars+=num_narrow_bars1; |
| 2310 | } |
| 2311 | |
| 2312 | // two short wide bars (4 wide in upgrade) |
| 2313 | Get("//composition[@name='forwardTOF_north']/mposY/@Y0",y0); |
| 2314 | Get("//composition[@name='forwardTOF_north']/mposY/@dY",dy); |
| 2315 | vector<double>tof_north; |
| 2316 | Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_north']/@X_Y_Z",tof_north); |
| 2317 | for (int k=num_bars;k<num_bars+num_single_end_bars1;k++){ |
| 2318 | y_tof.push_back(y0+tof_north[1]+dy*double(k-num_bars)); |
| 2319 | y_widths.push_back(dy); |
| 2320 | } |
| 2321 | num_bars+=num_single_end_bars1; |
| 2322 | |
| 2323 | // two narrower long bars - added by upgrade |
| 2324 | if(num_narrower_bars2 > 0) { |
| 2325 | Get("//composition[@name='forwardTOF_top3']/mposY/@Y0",y0); |
| 2326 | Get("//composition[@name='forwardTOF_top3']/mposY/@dY",dy); |
| 2327 | vector<double>tof_top3; |
| 2328 | Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_top3']/@X_Y_Z",tof_top3); |
| 2329 | for (int k=num_bars;k<num_bars+num_narrower_bars2;k++){ |
| 2330 | y_tof.push_back(y0+tof_top3[1]+dy*double(k-num_bars)); |
| 2331 | y_widths.push_back(dy); |
| 2332 | } |
| 2333 | num_bars+=num_narrow_bars2; |
| 2334 | } |
| 2335 | |
| 2336 | // two narrow long bars |
| 2337 | Get("//composition[@name='forwardTOF_top2']/mposY/@Y0",y0); |
| 2338 | Get("//composition[@name='forwardTOF_top2']/mposY/@dY",dy); |
| 2339 | vector<double>tof_top2; |
| 2340 | Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_top2']/@X_Y_Z",tof_top2); |
| 2341 | for (int k=num_bars;k<num_bars+num_narrow_bars2;k++){ |
| 2342 | y_tof.push_back(y0+tof_top2[1]+dy*double(k-num_bars)); |
| 2343 | y_widths.push_back(dy); |
| 2344 | } |
| 2345 | num_bars+=num_narrow_bars2; |
| 2346 | |
| 2347 | // Last 19 long bars |
| 2348 | Get("//composition[@name='forwardTOF_top1']/mposY/@Y0",y0); |
| 2349 | Get("//composition[@name='forwardTOF_top1']/mposY/@dY",dy); |
| 2350 | vector<double>tof_top1; |
| 2351 | Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_top1']/@X_Y_Z",tof_top1); |
| 2352 | for (int k=num_bars;k<num_bars+num_bars2;k++){ |
| 2353 | y_tof.push_back(y0+tof_top1[1]+dy*double(k-num_bars)); |
| 2354 | y_widths.push_back(dy); |
| 2355 | } |
| 2356 | num_bars+=num_bars2; |
Value stored to 'num_bars' is never read | |
| 2357 | |
| 2358 | /* |
| 2359 | // two more short wide bars - IGNORE FOR NOW, ASSUME SAME Y AS OTHER SINGLE ENDED |
| 2360 | Get("//composition[@name='forwardTOF_south']/mposY/@Y0",y0); |
| 2361 | Get("//composition[@name='forwardTOF_south']/mposY/@dY",dy); |
| 2362 | vector<double>tof_south; |
| 2363 | Get("//composition[@name='forwardTOF']/posXYZ[@volume='forwardTOF_south']/@X_Y_Z",tof_south); |
| 2364 | for (unsigned int k=45;k<47;k++){ |
| 2365 | y_tof.push_back(y0+tof_south[1]+dy*double(k-45)); |
| 2366 | } |
| 2367 | */ |
| 2368 | |
| 2369 | return true; |
| 2370 | } |
| 2371 | |
| 2372 | //--------------------------------- |
| 2373 | // GetTargetZ |
| 2374 | //--------------------------------- |
| 2375 | bool DGeometry::GetTargetZ(double &z_target) const |
| 2376 | { |
| 2377 | // Default to nominal center of GlueX target |
| 2378 | z_target=65.; |
| 2379 | |
| 2380 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 2381 | |
| 2382 | // Check GlueX target is defined |
| 2383 | bool gluex_target_exists = true; |
| 2384 | vector<double> xyz_vessel; |
| 2385 | vector<double> xyz_target; |
| 2386 | vector<double> xyz_detector; |
| 2387 | if(gluex_target_exists) gluex_target_exists = Get("//composition[@name='targetVessel']/posXYZ[@volume='targetTube']/@X_Y_Z", xyz_vessel); |
| 2388 | if(gluex_target_exists) gluex_target_exists = Get("//composition[@name='Target']/posXYZ[@volume='targetVessel']/@X_Y_Z", xyz_target); |
| 2389 | if(gluex_target_exists) gluex_target_exists = Get("//posXYZ[@volume='Target']/@X_Y_Z", xyz_detector); |
| 2390 | if(gluex_target_exists) { |
| 2391 | z_target = xyz_vessel[2] + xyz_target[2] + xyz_detector[2]; |
| 2392 | jgeom->SetVerbose(1); // reenable error messages |
| 2393 | return true; |
| 2394 | } |
| 2395 | |
| 2396 | // Check if CPP target is defined |
| 2397 | bool cpp_target_exists = true; |
| 2398 | vector<double> xyz_TGT0; |
| 2399 | vector<double> xyz_TARG; |
| 2400 | vector<double> xyz_TargetCPP; |
| 2401 | if(cpp_target_exists) cpp_target_exists = Get("//composition/posXYZ[@volume='TGT0']/@X_Y_Z", xyz_TGT0); |
| 2402 | if(cpp_target_exists) cpp_target_exists = Get("//composition/posXYZ[@volume='TARG']/@X_Y_Z", xyz_TARG); |
| 2403 | if(cpp_target_exists) cpp_target_exists = Get("//composition/posXYZ[@volume='TargetCPP']/@X_Y_Z", xyz_TargetCPP); |
| 2404 | if(cpp_target_exists) { |
| 2405 | z_target = xyz_TGT0[2] + xyz_TARG[2] + xyz_TargetCPP[2]; |
| 2406 | jgeom->SetVerbose(1); // reenable error messages |
| 2407 | return true; |
| 2408 | } |
| 2409 | |
| 2410 | // Check if PrimEx Be target is defined |
| 2411 | bool primex_target_exists = true; |
| 2412 | vector<double> xyz_BETG; |
| 2413 | if(primex_target_exists) primex_target_exists = Get("//composition[@name='targetVessel']/posXYZ[@volume='BETG']/@X_Y_Z", xyz_BETG); |
| 2414 | if(primex_target_exists) primex_target_exists = Get("//composition[@name='Target']/posXYZ[@volume='targetVessel']/@X_Y_Z", xyz_target); |
| 2415 | if(primex_target_exists) primex_target_exists = Get("//posXYZ[@volume='Target']/@X_Y_Z", xyz_detector); |
| 2416 | if(primex_target_exists) { |
| 2417 | |
| 2418 | z_target = xyz_BETG[2] + xyz_target[2] + xyz_detector[2]; |
| 2419 | |
| 2420 | //cout << " PrimEx Be targer selected. Z target = = " << z_target << endl; |
| 2421 | |
| 2422 | jgeom->SetVerbose(1); // reenable error messages |
| 2423 | return true; |
| 2424 | } |
| 2425 | |
| 2426 | // Check if SRC carbon foils are defined. Return the center of the ladder |
| 2427 | vector<double>ladder_xyz; |
| 2428 | if (Get("//composition[@name='targetVessel']/posXYZ[@volume='carbonLadder']/@X_Y_Z",ladder_xyz)){ |
| 2429 | Get("//composition[@name='Target']/posXYZ[@volume='targetVessel']/@X_Y_Z", xyz_target); |
| 2430 | Get("//posXYZ[@volume='Target']/@X_Y_Z", xyz_detector); |
| 2431 | |
| 2432 | vector<double>foil_pos; |
| 2433 | Get("//composition[@name='carbonLadder']/posXYZ[@volume='carbonTarget1']/@X_Y_Z",foil_pos); |
| 2434 | z_target=0.5*foil_pos[2]; |
| 2435 | Get("//composition[@name='carbonLadder']/posXYZ[@volume='carbonTarget8']/@X_Y_Z",foil_pos); |
| 2436 | z_target+=0.5*foil_pos[2]; |
| 2437 | z_target+=xyz_target[2]+xyz_detector[2]+ladder_xyz[2]; |
| 2438 | |
| 2439 | jgeom->SetVerbose(1); // reenable error messages |
| 2440 | return true; |
| 2441 | } |
| 2442 | |
| 2443 | // Only print warning for one thread whenever run number change |
| 2444 | static pthread_mutex_t empty_target_mutex = PTHREAD_MUTEX_INITIALIZER{ { 0, 0, 0, 0, 0, 0, 0, { 0, 0 } } }; |
| 2445 | static set<int> empty_target_runs_announced; |
| 2446 | |
| 2447 | // keep track of which runs we print out warnings for |
| 2448 | pthread_mutex_lock(&empty_target_mutex); |
| 2449 | bool empty_target_warning = false; |
| 2450 | if(empty_target_runs_announced.find(runnumber) == empty_target_runs_announced.end()){ |
| 2451 | empty_target_warning = true; |
| 2452 | empty_target_runs_announced.insert(runnumber); |
| 2453 | } |
| 2454 | pthread_mutex_unlock(&empty_target_mutex); |
| 2455 | |
| 2456 | if (empty_target_warning) |
| 2457 | 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 " << |
| 2458 | z_target << " cm" << endl; |
| 2459 | |
| 2460 | jgeom->SetVerbose(1); // reenable error messages |
| 2461 | |
| 2462 | return false; |
| 2463 | } |
| 2464 | |
| 2465 | //--------------------------------- |
| 2466 | // GetTargetLength |
| 2467 | //--------------------------------- |
| 2468 | bool DGeometry::GetTargetLength(double &target_length) const |
| 2469 | { |
| 2470 | |
| 2471 | target_length=0.; |
| 2472 | vector<double> zLength; |
| 2473 | jgeom->SetVerbose(0); // don't print error messages for optional detector elements |
| 2474 | |
| 2475 | // Regular GlueX target |
| 2476 | if (Get("//section[@name='Target']/pcon[@name='LIH2']/real[@name='length']/[@value]", zLength)){ |
| 2477 | target_length=zLength[0]; |
| 2478 | jgeom->SetVerbose(1); // reenable error messages |
| 2479 | return true; |
| 2480 | } |
| 2481 | |
| 2482 | // Liquid helium target |
| 2483 | if (Get("//section[@name='Target']/pcon[@name='LIHE']/real[@name='length']/[@value]", zLength)){ |
| 2484 | target_length=zLength[0]; |
| 2485 | jgeom->SetVerbose(1); // reenable error messages |
| 2486 | return true; |
| 2487 | } |
| 2488 | |
| 2489 | // Beryllium target |
| 2490 | vector<double>rio_z; |
| 2491 | if (Get("//section[@name='Target']/tubs[@name='BETG']/@Rio_Z",rio_z)){ |
| 2492 | target_length=rio_z[2]; |
| 2493 | jgeom->SetVerbose(1); // reenable error messages |
| 2494 | return true; |
| 2495 | } |
| 2496 | |
| 2497 | // Carbon foils for SRC; length is distance between first and last foil |
| 2498 | vector<double>foil_pos; |
| 2499 | if (Get("//composition[@name='carbonLadder']/posXYZ[@volume='carbonTarget1']/@X_Y_Z",foil_pos)){ |
| 2500 | target_length-=foil_pos[2]; |
| 2501 | Get("//composition[@name='carbonLadder']/posXYZ[@volume='carbonTarget8']/@X_Y_Z",foil_pos); |
| 2502 | target_length+=foil_pos[2]; |
| 2503 | jgeom->SetVerbose(1); // reenable error messages |
| 2504 | return true; |
| 2505 | } |
| 2506 | |
| 2507 | // Lead target for CPP |
| 2508 | if (Get("//section[@name='TargetCPP']/tubs[@name='TGT0']/@Rio_Z",rio_z)){ |
| 2509 | target_length=rio_z[2]; |
| 2510 | jgeom->SetVerbose(1); // reenable error messages |
| 2511 | return true; |
| 2512 | } |
| 2513 | |
| 2514 | jgeom->SetVerbose(1); // reenable error messages |
| 2515 | return false; |
| 2516 | } |
| 2517 | |
| 2518 | // Get vectors of positions and norm vectors for start counter from XML |
| 2519 | bool DGeometry::GetStartCounterGeom(vector<vector<DVector3> >&pos, |
| 2520 | vector<vector<DVector3> >&norm |
| 2521 | ) const |
| 2522 | { |
| 2523 | |
| 2524 | JCalibration *jcalib = dapp->GetJCalibration(runnumber); |
| 2525 | |
| 2526 | // Check if Start Counter geometry is present |
| 2527 | vector<double> sc_origin; |
| 2528 | bool got_sc = Get("//posXYZ[@volume='StartCntr']/@X_Y_Z", sc_origin); |
| 2529 | if(got_sc){ |
| 2530 | // Offset from beam center |
| 2531 | vector<double>sc_origin_delta; |
| 2532 | Get("//posXYZ[@volume='startCntr']/@X_Y_Z", sc_origin_delta); |
| 2533 | double dx=sc_origin_delta[0]; |
| 2534 | double dy=sc_origin_delta[1]; |
| 2535 | |
| 2536 | // z-position at upstream face of scintillators. |
| 2537 | double z0=sc_origin[2]; |
| 2538 | |
| 2539 | // Get rotation angles |
| 2540 | vector<double>sc_rot_angles; |
| 2541 | Get("//posXYZ[@volume='startCntr']/@rot", sc_rot_angles); |
| 2542 | double ThetaX=sc_rot_angles[0]*M_PI3.14159265358979323846/180.; |
| 2543 | double ThetaY=sc_rot_angles[1]*M_PI3.14159265358979323846/180.; |
| 2544 | Get("//posXYZ[@volume='StartCntr']/@rot", sc_rot_angles); |
| 2545 | //double ThetaX=sc_rot_angles[0]*M_PI/180.; |
| 2546 | //double ThetaY=sc_rot_angles[1]*M_PI/180.; |
| 2547 | double ThetaZ=sc_rot_angles[2]*M_PI3.14159265358979323846/180.; |
| 2548 | |
| 2549 | // Get overall alignment shifts from CCDB |
| 2550 | map<string,double> sc_global_offsets; |
| 2551 | if (jcalib->Get("START_COUNTER/global_alignment_parms",sc_global_offsets)==false) { |
| 2552 | // translations |
| 2553 | dx += sc_global_offsets["SC_ALIGN_X"]; |
| 2554 | dy += sc_global_offsets["SC_ALIGN_Y"]; |
| 2555 | z0 += sc_global_offsets["SC_ALIGN_Z"]; |
| 2556 | |
| 2557 | // rotations |
| 2558 | ThetaX += sc_global_offsets["SC_ALIGN_ROTX"]*M_PI3.14159265358979323846/180.; |
| 2559 | ThetaY += sc_global_offsets["SC_ALIGN_ROTY"]*M_PI3.14159265358979323846/180.; |
| 2560 | ThetaZ += sc_global_offsets["SC_ALIGN_ROTZ"]*M_PI3.14159265358979323846/180.; |
| 2561 | } |
| 2562 | |
| 2563 | double num_paddles; |
| 2564 | Get("//mposPhi[@volume='STRC']/@ncopy",num_paddles); |
| 2565 | double dSCdphi = M_TWO_PI6.28318530717958647692/num_paddles; |
| 2566 | |
| 2567 | vector<vector<double> > sc_rioz; |
| 2568 | GetMultiple("//pgon[@name='STRC']/polyplane/@Rio_Z", sc_rioz); |
| 2569 | |
| 2570 | // Get individual paddle alignment parameters |
| 2571 | vector< map<string,double> > sc_paddle_offsets; |
| 2572 | bool loaded_paddle_offsets = false; |
| 2573 | if (jcalib->Get("START_COUNTER/paddle_alignment_parms",sc_paddle_offsets)==false) |
| 2574 | loaded_paddle_offsets = true; |
| 2575 | |
| 2576 | // Create vectors of positions and normal vectors for each paddle |
| 2577 | for (unsigned int i=0;i<30;i++){ |
| 2578 | double phi=ThetaZ+dSCdphi*(double(i)+0.5); |
| 2579 | double sinphi=sin(phi); |
| 2580 | double cosphi=cos(phi); |
| 2581 | double r=0.5*(sc_rioz[0][0]+sc_rioz[0][1]); |
| 2582 | DVector3 oldray; |
| 2583 | // Rotate by phi and take into account the tilt |
| 2584 | double x=r*cosphi;//+dx; |
| 2585 | double y=r*sinphi;//+dy; |
| 2586 | DVector3 ray(x,y,sc_rioz[0][2]); |
| 2587 | ray.RotateX(ThetaX); |
| 2588 | ray.RotateY(ThetaY); |
| 2589 | |
| 2590 | // Create stl-vectors to store positions and norm vectors |
| 2591 | vector<DVector3>posvec; |
| 2592 | vector<DVector3>dirvec; |
| 2593 | // Loop over radial/z positions describing start counter geometry from xml |
| 2594 | for(unsigned int k = 1; k < sc_rioz.size(); ++k){ |
| 2595 | oldray=ray; |
| 2596 | r=0.5*(sc_rioz[k][0]+sc_rioz[k][1]); |
| 2597 | // Point in plane of scintillator |
| 2598 | x=r*cosphi;//+dx; |
| 2599 | y=r*sinphi;//+dy; |
| 2600 | ray.SetXYZ(x,y,sc_rioz[k][2]); |
| 2601 | ray.RotateX(ThetaX); |
| 2602 | ray.RotateY(ThetaY); |
| 2603 | // Apply alignment parameters |
| 2604 | if(loaded_paddle_offsets) { |
| 2605 | // allow for a maximum extent of the paddle in z |
| 2606 | double max_z = sc_paddle_offsets[i]["SC_MAX_Z"]; |
| 2607 | if(ray.Z() > max_z) { |
| 2608 | ray.SetZ(max_z); |
| 2609 | } |
| 2610 | // allow for a modification of the bend angle of the paddle (18.5 deg from horizontal) |
| 2611 | // this should just be a perturbation around this angle, so assume a linear interpolation |
| 2612 | double delta_theta = sc_paddle_offsets[i]["SC_CURVE_THETA"]; // in degrees, r of curvature = 120 cm |
| 2613 | ray.SetX(ray.X()+delta_theta*1.65); // 1 degree ~ 1.65 cm in x |
| 2614 | ray.SetY(ray.Y()-delta_theta*0.55); // 1 degree ~ 0.55 cm in y |
| 2615 | } |
| 2616 | // Second point in the plane of the scintillator |
| 2617 | DVector3 ray2(r,10.,sc_rioz[k][2]); |
| 2618 | ray2.RotateZ(phi+0.5*dSCdphi*(1.+1./15.*((i>14)?29-i:i))); |
| 2619 | ray2.RotateX(ThetaX); |
| 2620 | ray2.RotateY(ThetaY); |
| 2621 | // Compute normal vector to plane |
| 2622 | DVector3 dir=(ray-oldray).Cross(ray2-oldray); |
| 2623 | dir.SetMag(1.); |
| 2624 | dirvec.push_back(dir); |
| 2625 | posvec.push_back(DVector3(oldray.X()+dx,oldray.Y()+dy,oldray.Z()+z0)); |
| 2626 | } |
| 2627 | posvec.push_back(DVector3(ray.X()+dx,ray.Y()+dy,ray.Z()+z0)); //SAVE THE ENDPOINT OF THE LAST PLANE |
| 2628 | pos.push_back(posvec); |
| 2629 | norm.push_back(dirvec); |
| 2630 | |
| 2631 | posvec.clear(); |
| 2632 | dirvec.clear(); |
| 2633 | } |
| 2634 | |
| 2635 | } |
| 2636 | return got_sc; |
| 2637 | } |
| 2638 |