GeomBuilder/src/MagGeoBuilderFromDDD.cc

0001 /*
0002  *  See header file for a description of this class.
0003  *
0004  *  \author N. Amapane - INFN Torino
0005  */
0006
0007 #include "MagGeoBuilderFromDDD.h"
0008 #include "volumeHandle.h"
0009 #include "bSlab.h"
0010 #include "bRod.h"
0011 #include "bSector.h"
0012 #include "bLayer.h"
0013 #include "eSector.h"
0014 #include "eLayer.h"
0015 #include "FakeInterpolator.h"
0016
0017 #include "MagneticField/Layers/interface/MagBLayer.h"
0018 #include "MagneticField/Layers/interface/MagESector.h"
0019
0020 #include "FWCore/ParameterSet/interface/FileInPath.h"
0021
0022 #include "DetectorDescription/Core/interface/DDCompactView.h"
0023 #include "DetectorDescription/Core/interface/DDFilteredView.h"
0024 #include "DetectorDescription/Core/interface/DDFilter.h"
0025
0026 #include "Utilities/BinningTools/interface/ClusterizingHistogram.h"
0027
0028 #include "MagneticField/Interpolation/interface/MagProviderInterpol.h"
0029 #include "MagneticField/Interpolation/interface/MFGridFactory.h"
0030 #include "MagneticField/Interpolation/interface/MFGrid.h"
0031
0032 #include "MagneticField/VolumeGeometry/interface/MagVolume6Faces.h"
0033 #include "MagneticField/VolumeGeometry/interface/MagExceptions.h"
0034
0035 #include "DataFormats/GeometryVector/interface/Pi.h"
0036 #include "DataFormats/Math/interface/GeantUnits.h"
0037
0038 #include "FWCore/Utilities/interface/Exception.h"
0039 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0040
0041 #include <string>
0042 #include <vector>
0043 #include <iostream>
0044 #include <sstream>
0045 #include <algorithm>
0046 #include <iterator>
0047 #include <map>
0048 #include <set>
0049 #include <iomanip>
0050 #include <boost/algorithm/string/replace.hpp>
0051 #include "Utilities/General/interface/precomputed_value_sort.h"
0052
0053 using namespace std;
0054 using namespace magneticfield;
0055
0056 MagGeoBuilderFromDDD::MagGeoBuilderFromDDD(string tableSet_, int geometryVersion_, bool debug_)
0057     : tableSet(tableSet_), geometryVersion(geometryVersion_), theGridFiles(nullptr), debug(debug_) {
0058   if (debug)
0059     cout << "Constructing a MagGeoBuilderFromDDD" << endl;
0060 }
0061
0062 MagGeoBuilderFromDDD::~MagGeoBuilderFromDDD() {
0063   for (handles::const_iterator i = bVolumes.begin(); i != bVolumes.end(); ++i) {
0064     delete (*i);
0065   }
0066
0067   for (handles::const_iterator i = eVolumes.begin(); i != eVolumes.end(); ++i) {
0068     delete (*i);
0069   }
0070 }
0071
0072 void MagGeoBuilderFromDDD::summary(handles& volumes) {
0073   // The final countdown.
0074   int ivolumes = volumes.size();  // number of volumes
0075   int isurfaces = ivolumes * 6;   // number of individual surfaces
0076   int iassigned = 0;              // How many have been assigned
0077   int iunique = 0;                // number of unique surfaces
0078   int iref_ass = 0;
0079   int iref_nass = 0;
0080
0081   set<const void*> ptrs;
0082
0083   handles::const_iterator first = volumes.begin();
0084   handles::const_iterator last = volumes.end();
0085
0086   for (handles::const_iterator i = first; i != last; ++i) {
0087     if (int((*i)->shape()) > 4)
0088       continue;  // FIXME: implement test for missing shapes...
0089     for (int side = 0; side < 6; ++side) {
0090       int references = (*i)->references(side);
0091       if ((*i)->isPlaneMatched(side)) {
0092         ++iassigned;
0093         bool firstOcc = (ptrs.insert(&((*i)->surface(side)))).second;
0094         if (firstOcc)
0095           iref_ass += references;
0096         if (references < 2) {
0097           cout << "*** Only 1 ref, vol: " << (*i)->volumeno << " # " << (*i)->copyno << " side: " << side << endl;
0098         }
0099       } else {
0100         iref_nass += references;
0101         if (references > 1) {
0102           cout << "*** Ref_nass >1 " << endl;
0103         }
0104       }
0105     }
0106   }
0107   iunique = ptrs.size();
0108
0109   cout << "    volumes   " << ivolumes << endl
0110        << "    surfaces  " << isurfaces << endl
0111        << "    assigned  " << iassigned << endl
0112        << "    unique    " << iunique << endl
0113        << "    iref_ass  " << iref_ass << endl
0114        << "    iref_nass " << iref_nass << endl;
0115 }
0116
0117 void MagGeoBuilderFromDDD::build(const DDCompactView& cpva) {
0118   //    DDCompactView cpv;
0119   DDExpandedView fv(cpva);
0120
0121   if (debug)
0122     cout << "**********************************************************" << endl;
0123
0124   // The actual field interpolators
0125   map<string, MagProviderInterpol*> bInterpolators;
0126   map<string, MagProviderInterpol*> eInterpolators;
0127
0128   // Counter of different volumes
0129   int bVolCount = 0;
0130   int eVolCount = 0;
0131
0132   if (fv.logicalPart().name().name() != "MAGF") {
0133     std::string topNodeName(fv.logicalPart().name().name());
0134
0135     //see if one of the children is MAGF
0136     bool doSubDets = fv.firstChild();
0137
0138     bool go = true;
0139     while (go && doSubDets) {
0140       if (fv.logicalPart().name().name() == "MAGF")
0141         break;
0142       else
0143         go = fv.nextSibling();
0144     }
0145     if (!go) {
0146       throw cms::Exception("NoMAGFinDDD")
0147           << " Neither the top node, nor any child node of the DDCompactView is \"MAGF\" but the top node is instead \""
0148           << topNodeName << "\"";
0149     }
0150   }
0151   // Loop over MAGF volumes and create volumeHandles.
0152   if (debug) {
0153     cout << endl
0154          << "*** MAGF: " << fv.geoHistory() << endl
0155          << "translation: " << fv.translation() << endl
0156          << " rotation: " << fv.rotation() << endl;
0157   }
0158
0159   bool doSubDets = fv.firstChild();
0160   while (doSubDets) {
0161     string name = fv.logicalPart().name().name();
0162     if (debug)
0163       cout << endl << "Name: " << name << endl << "      " << fv.geoHistory() << endl;
0164
0165     // FIXME: single-volyme cylinders - this is feature has been removed
0166     //        and should be revisited.
0167     //    bool mergeCylinders=false;
0168
0169     // If true, In the barrel, cylinders sectors will be skipped to build full
0170     // cylinders out of sector copyno #1.
0171     bool expand = false;
0172
0173     //     if (mergeCylinders) {
0174     //       if (name == "V_ZN_1"
0175     //    || name == "V_ZN_2") {
0176     //  if (debug && fv.logicalPart().solid().shape()!=ddtubs) {
0177     //    cout << "ERROR: MagGeoBuilderFromDDD::build: volume " << name
0178     //         << " should be a cylinder" << endl;
0179     //  }
0180     //  if(fv.copyno()==1) {
0181     //    expand = true;
0182     //  } else {
0183     //    //cout << "... to be skipped: "
0184     //    //     << name << " " << fv.copyno() << endl;
0185     //  }
0186     //       }
0187     //     }
0188
0189     volumeHandle* v = new volumeHandle(fv, expand, debug);
0190
0191     if (theGridFiles != nullptr) {
0192       int key = (v->volumeno) * 100 + v->copyno;
0193       TableFileMap::const_iterator itable = theGridFiles->find(key);
0194       if (itable == theGridFiles->end()) {
0195         key = (v->volumeno) * 100;
0196         itable = theGridFiles->find(key);
0197       }
0198
0199       if (itable != theGridFiles->end()) {
0200         string magFile = (*itable).second.first;
0201         stringstream conv;
0202         string svol, ssec;
0203         conv << setfill('0') << setw(3) << v->volumeno << " " << setw(2)
0204              << v->copyno;  // volume assumed to have 0s padding to 3 digits; sector assumed to have 0s padding to 2 digits
0205         conv >> svol >> ssec;
0206         boost::replace_all(magFile, "[v]", svol);
0207         boost::replace_all(magFile, "[s]", ssec);
0208         int masterSector = (*itable).second.second;
0209         if (masterSector == 0)
0210           masterSector = v->copyno;
0211         v->magFile = magFile;
0212         v->masterSector = masterSector;
0213       } else {
0214         edm::LogError("MagGeoBuilderFromDDDbuild")
0215             << "ERROR: no table spec found for V " << v->volumeno << ":" << v->copyno;
0216       }
0217     }
0218
0219     // Select volumes, build volume handles.
0220     float Z = v->center().z();
0221     float R = v->center().perp();
0222
0223     // v 85l: Barrel is everything up to |Z| = 661.0, excluding
0224     // volume #7, centered at 6477.5
0225     // v 1103l: same numbers work fine. #16 instead of #7, same coords;
0226     // see comment below for V6,7
0227     //ASSUMPTION: no misalignment is applied to mag volumes.
0228     //FIXME: implement barrel/endcap flags as DDD SpecPars.
0229     if ((fabs(Z) < 647. || (R > 350. && fabs(Z) < 662.)) &&
0230         !(fabs(Z) > 480 && R < 172)  // in 1103l we place V_6 and V_7 in the
0231                                      // endcaps to preserve nice layer structure
0232                                      // in the barrel. This does not hurt in v85l
0233                                      // where there is a single V1
0234     ) {                              // Barrel
0235       if (debug)
0236         cout << " (Barrel)" << endl;
0237       bVolumes.push_back(v);
0238
0239       // Build the interpolator of the "master" volume (the one which is
0240       // not replicated in phi)
0241       // ASSUMPTION: copyno == sector.
0242       if (v->copyno == v->masterSector) {
0243         buildInterpolator(v, bInterpolators);
0244         ++bVolCount;
0245       }
0246     } else {  // Endcaps
0247       if (debug)
0248         cout << " (Endcaps)" << endl;
0249       eVolumes.push_back(v);
0250       if (v->copyno == v->masterSector) {
0251         buildInterpolator(v, eInterpolators);
0252         ++eVolCount;
0253       }
0254     }
0255
0256     doSubDets = fv.nextSibling();  // end of loop over MAGF
0257   }
0258
0259   if (debug) {
0260     cout << "Number of volumes (barrel): " << bVolumes.size() << endl
0261          << "Number of volumes (endcap): " << eVolumes.size() << endl;
0262     cout << "**********************************************************" << endl;
0263   }
0264
0265   // Now all volumeHandles are there, and parameters for each of the planes
0266   // are calculated.
0267
0268   //----------------------------------------------------------------------
0269   // Print summary information
0270
0271   if (debug) {
0272     cout << "-----------------------" << endl;
0273     cout << "SUMMARY: Barrel " << endl;
0274     summary(bVolumes);
0275
0276     cout << endl << "SUMMARY: Endcaps " << endl;
0277     summary(eVolumes);
0278     cout << "-----------------------" << endl;
0279   }
0280
0281   //----------------------------------------------------------------------
0282   // Find barrel layers.
0283
0284   vector<bLayer> layers;  // the barrel layers
0285   precomputed_value_sort(bVolumes.begin(), bVolumes.end(), ExtractRN());
0286
0287   // Find the layers (in R)
0288   const float resolution = 1.;  // cm
0289   float rmin = bVolumes.front()->RN() - resolution;
0290   float rmax = bVolumes.back()->RN() + resolution;
0291   ClusterizingHistogram hisR(int((rmax - rmin) / resolution) + 1, rmin, rmax);
0292
0293   if (debug)
0294     cout << " R layers: " << rmin << " " << rmax << endl;
0295
0296   handles::const_iterator first = bVolumes.begin();
0297   handles::const_iterator last = bVolumes.end();
0298
0299   for (handles::const_iterator i = first; i != last; ++i) {
0300     hisR.fill((*i)->RN());
0301   }
0302   vector<float> rClust = hisR.clusterize(resolution);
0303
0304   handles::const_iterator ringStart = first;
0305   handles::const_iterator separ = first;
0306
0307   for (unsigned int i = 0; i < rClust.size() - 1; ++i) {
0308     if (debug)
0309       cout << " Layer at RN = " << rClust[i];
0310     float rSepar = (rClust[i] + rClust[i + 1]) / 2.f;
0311     while ((*separ)->RN() < rSepar)
0312       ++separ;
0313
0314     bLayer thislayer(ringStart, separ, debug);
0315     layers.push_back(thislayer);
0316     ringStart = separ;
0317   }
0318   {
0319     if (debug)
0320       cout << " Layer at RN = " << rClust.back();
0321     bLayer thislayer(separ, last, debug);
0322     layers.push_back(thislayer);
0323   }
0324
0325   if (debug)
0326     cout << "Barrel: Found " << rClust.size() << " clusters in R, " << layers.size() << " layers " << endl << endl;
0327
0328   //----------------------------------------------------------------------
0329   // Find endcap sectors
0330   vector<eSector> sectors;  // the endcap sectors
0331
0332   // Find the number of sectors (should be 12 or 24 depending on the geometry model)
0333   float phireso = 0.05;  // rad
0334   ClusterizingHistogram hisPhi(int((Geom::ftwoPi()) / phireso) + 1, -Geom::fpi(), Geom::fpi());
0335
0336   for (handles::const_iterator i = eVolumes.begin(); i != eVolumes.end(); ++i) {
0337     hisPhi.fill((*i)->minPhi());
0338   }
0339   vector<float> phiClust = hisPhi.clusterize(phireso);
0340   int nESectors = phiClust.size();
0341   if (debug && (nESectors % 12) != 0)
0342     cout << "ERROR: unexpected # of sectors: " << nESectors << endl;
0343
0344   //Sort in phi
0345   precomputed_value_sort(eVolumes.begin(), eVolumes.end(), ExtractPhi());
0346
0347   // Handle the -pi/pi boundary: volumes crossing it could be half at the begin and half at end of the sorted list.
0348   // So, check if any of the volumes that should belong to the first bin (at -phi) are at the end of the list:
0349   float lastBinPhi = phiClust.back();
0350   handles::reverse_iterator ri = eVolumes.rbegin();
0351   while ((*ri)->center().phi() > lastBinPhi) {
0352     ++ri;
0353   }
0354   if (ri != eVolumes.rbegin()) {
0355     // ri points to the first element that is within the last bin.
0356     // We need to move the following element (ie ri.base()) to the beginning of the list,
0357     handles::iterator newbeg = ri.base();
0358     rotate(eVolumes.begin(), newbeg, eVolumes.end());
0359   }
0360
0361   //Group volumes in sectors
0362   int offset = eVolumes.size() / nESectors;
0363   for (int i = 0; i < nESectors; ++i) {
0364     if (debug) {
0365       cout << " Sector at phi = " << (*(eVolumes.begin() + ((i)*offset)))->center().phi() << endl;
0366       // Additional x-check: sectors are expected to be made by volumes with the same copyno
0367       int secCopyNo = -1;
0368       for (handles::const_iterator iv = eVolumes.begin() + ((i)*offset); iv != eVolumes.begin() + ((i + 1) * offset);
0369            ++iv) {
0370         if (secCopyNo >= 0 && (*iv)->copyno != secCopyNo)
0371           cout << "ERROR: volume copyno" << (*iv)->name << ":" << (*iv)->copyno
0372                << " differs from others in same sectors " << secCopyNo << endl;
0373         secCopyNo = (*iv)->copyno;
0374       }
0375     }
0376
0377     sectors.push_back(eSector(eVolumes.begin() + ((i)*offset), eVolumes.begin() + ((i + 1) * offset), debug));
0378   }
0379
0380   if (debug)
0381     cout << "Endcap: Found " << sectors.size() << " sectors " << endl;
0382
0383   //----------------------------------------------------------------------
0384   // Match surfaces.
0385
0386   //  cout << "------------------" << endl << "Now associating planes..." << endl;
0387
0388   //   // Loop on layers
0389   //   for (vector<bLayer>::const_iterator ilay = layers.begin();
0390   //        ilay!= layers.end(); ++ilay) {
0391   //     cout << "On Layer: " << ilay-layers.begin() << " RN: " << (*ilay).RN()
0392   //     <<endl;
0393
0394   //     // Loop on wheels
0395   //     for (vector<bWheel>::const_iterator iwheel = (*ilay).wheels.begin();
0396   //     iwheel != (*ilay).wheels.end(); ++iwheel) {
0397   //       cout << "  On Wheel: " << iwheel- (*ilay).wheels.begin()<< " Z: "
0398   //       << (*iwheel).minZ() << " " << (*iwheel).maxZ() << " "
0399   //       << ((*iwheel).minZ()+(*iwheel).maxZ())/2. <<endl;
0400
0401   //       // Loop on sectors.
0402   //       for (int isector = 0; isector<12; ++isector) {
0403   //    // FIXME: create new constructor...
0404   //    bSectorNavigator navy(layers,
0405   //                  ilay-layers.begin(),
0406   //                  iwheel-(*ilay).wheels.begin(),isector);
0407
0408   //    const bSector & isect = (*iwheel).sector(isector);
0409
0410   //    isect.matchPlanes(navy); //FIXME refcount
0411   //       }
0412   //     }
0413   //   }
0414
0415   //----------------------------------------------------------------------
0416   // Build MagVolumes and the MagGeometry hierarchy.
0417
0418   //--- Barrel
0419
0420   // Build MagVolumes and associate interpolators to them
0421   buildMagVolumes(bVolumes, bInterpolators);
0422
0423   // Build MagBLayers
0424   for (vector<bLayer>::const_iterator ilay = layers.begin(); ilay != layers.end(); ++ilay) {
0425     mBLayers.push_back((*ilay).buildMagBLayer());
0426   }
0427
0428   if (debug) {
0429     cout << "*** BARREL ********************************************" << endl
0430          << "Number of different volumes   = " << bVolCount << endl
0431          << "Number of interpolators built = " << bInterpolators.size() << endl
0432          << "Number of MagBLayers built    = " << mBLayers.size() << endl;
0433
0434     testInside(bVolumes);  // FIXME: all volumes should be checked in one go.
0435   }
0436
0437   //--- Endcap
0438   // Build MagVolumes  and associate interpolators to them
0439   buildMagVolumes(eVolumes, eInterpolators);
0440
0441   // Build the MagESectors
0442   for (vector<eSector>::const_iterator isec = sectors.begin(); isec != sectors.end(); ++isec) {
0443     mESectors.push_back((*isec).buildMagESector());
0444   }
0445
0446   if (debug) {
0447     cout << "*** ENDCAP ********************************************" << endl
0448          << "Number of different volumes   = " << eVolCount << endl
0449          << "Number of interpolators built = " << eInterpolators.size() << endl
0450          << "Number of MagESector built    = " << mESectors.size() << endl;
0451
0452     testInside(eVolumes);  // FIXME: all volumes should be checked in one go.
0453   }
0454 }
0455
0456 void MagGeoBuilderFromDDD::buildMagVolumes(const handles& volumes, map<string, MagProviderInterpol*>& interpolators) {
0457   // Build all MagVolumes setting the MagProviderInterpol
0458   for (handles::const_iterator vol = volumes.begin(); vol != volumes.end(); ++vol) {
0459     const MagProviderInterpol* mp = nullptr;
0460     if (interpolators.find((*vol)->magFile) != interpolators.end()) {
0461       mp = interpolators[(*vol)->magFile];
0462     } else {
0463       edm::LogError("MagGeoBuilder") << "No interpolator found for file " << (*vol)->magFile
0464                                      << " vol: " << (*vol)->volumeno << "\n"
0465                                      << interpolators.size() << endl;
0466     }
0467
0468     // Search for [volume,sector] in the list of scaling factors; sector = 0 handled as wildcard
0469     // ASSUMPTION: copyno == sector.
0470     int key = ((*vol)->volumeno) * 100 + (*vol)->copyno;
0471     map<int, double>::const_iterator isf = theScalingFactors.find(key);
0472     if (isf == theScalingFactors.end()) {
0473       key = ((*vol)->volumeno) * 100;
0474       isf = theScalingFactors.find(key);
0475     }
0476
0477     double sf = 1.;
0478     if (isf != theScalingFactors.end()) {
0479       sf = (*isf).second;
0480
0481       edm::LogInfo("MagGeoBuilder") << "Applying scaling factor " << sf << " to " << (*vol)->volumeno << "["
0482                                     << (*vol)->copyno << "] (key:" << key << ")" << endl;
0483     }
0484
0485     const GloballyPositioned<float>* gpos = (*vol)->placement();
0486     (*vol)->magVolume = new MagVolume6Faces(gpos->position(), gpos->rotation(), (*vol)->sides(), mp, sf);
0487
0488     if ((*vol)->copyno == (*vol)->masterSector) {
0489       (*vol)->magVolume->ownsFieldProvider(true);
0490     }
0491
0492     (*vol)->magVolume->setIsIron((*vol)->isIron());
0493
0494     // The name and sector of the volume are saved for debug purposes only. They may be removed at some point...
0495     (*vol)->magVolume->volumeNo = (*vol)->volumeno;
0496     (*vol)->magVolume->copyno = (*vol)->copyno;
0497   }
0498 }
0499
0500 void MagGeoBuilderFromDDD::buildInterpolator(const volumeHandle* vol,
0501                                              map<string, MagProviderInterpol*>& interpolators) {
0502   // Phi of the master sector
0503   double masterSectorPhi = (vol->masterSector - 1) * Geom::pi() / 6.;
0504
0505   if (debug) {
0506     cout << "Building interpolator from " << vol->volumeno << " copyno " << vol->copyno << " at " << vol->center()
0507          << " phi: " << vol->center().phi() << " file: " << vol->magFile << " master : " << vol->masterSector << endl;
0508
0509     if (fabs(vol->center().phi() - masterSectorPhi) > Geom::pi() / 9.) {
0510       cout << "***WARNING wrong sector? " << endl;
0511     }
0512   }
0513
0514   if (tableSet == "fake" || vol->magFile == "fake") {
0515     interpolators[vol->magFile] = new magneticfield::FakeInterpolator();
0516     return;
0517   }
0518
0519   string fullPath;
0520
0521   try {
0522     edm::FileInPath mydata("MagneticField/Interpolation/data/" + tableSet + "/" + vol->magFile);
0523     fullPath = mydata.fullPath();
0524   } catch (edm::Exception& exc) {
0525     cerr << "MagGeoBuilderFromDDD: exception in reading table; " << exc.what() << endl;
0526     if (!debug)
0527       throw;
0528     return;
0529   }
0530
0531   try {
0532     if (vol->toExpand()) {
0533       //FIXME: see discussion on mergeCylinders above.
0534       //       interpolators[vol->magFile] =
0535       //    MFGridFactory::build( fullPath, *(vol->placement()), vol->minPhi(), vol->maxPhi());
0536     } else {
0537       // If the table is in "local" coordinates, must create a reference
0538       // frame that is appropriately rotated along the CMS Z axis.
0539
0540       GloballyPositioned<float> rf = *(vol->placement());
0541
0542       if (vol->masterSector != 1) {
0543         typedef Basic3DVector<float> Vector;
0544
0545         GloballyPositioned<float>::RotationType rot(Vector(0, 0, 1), -masterSectorPhi);
0546         Vector vpos(vol->placement()->position());
0547
0548         rf = GloballyPositioned<float>(GloballyPositioned<float>::PositionType(rot.multiplyInverse(vpos)),
0549                                        vol->placement()->rotation() * rot);
0550       }
0551
0552       interpolators[vol->magFile] = MFGridFactory::build(fullPath, rf);
0553     }
0554   } catch (MagException& exc) {
0555     cout << exc.what() << endl;
0556     interpolators.erase(vol->magFile);
0557     if (!debug)
0558       throw;
0559     return;
0560   }
0561
0562   if (debug) {
0563     // Check that all grid points of the interpolator are inside the volume.
0564     const MagVolume6Faces tempVolume(
0565         vol->placement()->position(), vol->placement()->rotation(), vol->sides(), interpolators[vol->magFile]);
0566
0567     const MFGrid* grid = dynamic_cast<const MFGrid*>(interpolators[vol->magFile]);
0568     if (grid != nullptr) {
0569       Dimensions sizes = grid->dimensions();
0570       cout << "Grid has 3 dimensions "
0571            << " number of nodes is " << sizes.w << " " << sizes.h << " " << sizes.d << endl;
0572
0573       const double tolerance = 0.03;
0574
0575       size_t dumpCount = 0;
0576       for (int j = 0; j < sizes.h; j++) {
0577         for (int k = 0; k < sizes.d; k++) {
0578           for (int i = 0; i < sizes.w; i++) {
0579             MFGrid::LocalPoint lp = grid->nodePosition(i, j, k);
0580             if (!tempVolume.inside(lp, tolerance)) {
0581               if (++dumpCount < 2) {
0582                 MFGrid::GlobalPoint gp = tempVolume.toGlobal(lp);
0583                 cout << "GRID ERROR: " << i << " " << j << " " << k << " local: " << lp << " global: " << gp
0584                      << " R= " << gp.perp() << " phi=" << gp.phi() << endl;
0585               }
0586             }
0587           }
0588         }
0589       }
0590
0591       cout << "Volume:" << vol->volumeno << " : Number of grid points outside the MagVolume: " << dumpCount << "/"
0592            << sizes.w * sizes.h * sizes.d << endl;
0593     }
0594   }
0595 }
0596
0597 void MagGeoBuilderFromDDD::testInside(handles& volumes) {
0598   // test inside() for all volumes.
0599   cout << "--------------------------------------------------" << endl;
0600   cout << " inside(center) test" << endl;
0601   for (handles::const_iterator vol = volumes.begin(); vol != volumes.end(); ++vol) {
0602     for (handles::const_iterator i = volumes.begin(); i != volumes.end(); ++i) {
0603       if ((*i) == (*vol))
0604         continue;
0605       //if ((*i)->magVolume == 0) continue;
0606       if ((*i)->magVolume->inside((*vol)->center())) {
0607         cout << "*** ERROR: center of V " << (*vol)->volumeno << ":" << (*vol)->copyno << " is inside V "
0608              << (*i)->volumeno << ":" << (*i)->copyno << endl;
0609       }
0610     }
0611
0612     if ((*vol)->magVolume->inside((*vol)->center())) {
0613       cout << "V " << (*vol)->volumeno << " OK " << endl;
0614     } else {
0615       cout << "*** ERROR: center of volume is not inside it, " << (*vol)->volumeno << endl;
0616     }
0617   }
0618   cout << "--------------------------------------------------" << endl;
0619 }
0620
0621 vector<MagBLayer*> MagGeoBuilderFromDDD::barrelLayers() const { return mBLayers; }
0622
0623 vector<MagESector*> MagGeoBuilderFromDDD::endcapSectors() const { return mESectors; }
0624
0625 vector<MagVolume6Faces*> MagGeoBuilderFromDDD::barrelVolumes() const {
0626   vector<MagVolume6Faces*> v;
0627   v.reserve(bVolumes.size());
0628   for (handles::const_iterator i = bVolumes.begin(); i != bVolumes.end(); ++i) {
0629     v.push_back((*i)->magVolume);
0630   }
0631   return v;
0632 }
0633
0634 vector<MagVolume6Faces*> MagGeoBuilderFromDDD::endcapVolumes() const {
0635   vector<MagVolume6Faces*> v;
0636   v.reserve(eVolumes.size());
0637   for (handles::const_iterator i = eVolumes.begin(); i != eVolumes.end(); ++i) {
0638     v.push_back((*i)->magVolume);
0639   }
0640   return v;
0641 }
0642
0643 float MagGeoBuilderFromDDD::maxR() const {
0644   //FIXME: should get it from the actual geometry
0645   return 900.;
0646 }
0647
0648 float MagGeoBuilderFromDDD::maxZ() const {
0649   //FIXME: should get it from the actual geometry
0650   if (geometryVersion >= 160812)
0651     return 2400.;
0652   else if (geometryVersion >= 120812)
0653     return 2000.;
0654   else
0655     return 1600.;
0656 }
0657
0658 void MagGeoBuilderFromDDD::setScaling(const std::vector<int>& keys, const std::vector<double>& values) {
0659   if (keys.size() != values.size()) {
0660     throw cms::Exception("InvalidParameter")
0661         << "Invalid field scaling parameters 'scalingVolumes' and 'scalingFactors' ";
0662   }
0663   for (unsigned int i = 0; i < keys.size(); ++i) {
0664     theScalingFactors[keys[i]] = values[i];
0665   }
0666 }
0667
0668 void MagGeoBuilderFromDDD::setGridFiles(const TableFileMap& gridFiles) { theGridFiles = &gridFiles; }