From 957d80d086fb0b5622f59510ed3b08429459e353 Mon Sep 17 00:00:00 2001 From: Carl Vuosalo Date: Wed, 3 Jul 2019 01:20:47 +0200 Subject: [PATCH] Condensed commit for mf geom builder migration to DD4hep --- .../DDCMS/data/cms-mf-geometry.xml | 1 + .../DDCMS/plugins/BuildFile.xml | 1 + .../plugins/DDCompactViewMFESProducer.cc | 70 ++ .../GeomBuilder/plugins/BuildFile.xml | 7 + .../plugins/dd4hep/MagGeoBuilder.cc | 631 ++++++++++++++++++ .../plugins/dd4hep/MagGeoBuilder.h | 98 +++ .../dd4hep/VolBasedMagFieldESProducerNewDD.cc | 114 ++++ .../plugins/dd4hep/volumeHandle.cc | 207 ++++++ .../GeomBuilder/plugins/dd4hep/volumeHandle.h | 60 ++ .../GeomBuilder/src/BaseVolumeHandle.cc | 221 ++++++ .../GeomBuilder/src/BaseVolumeHandle.h | 219 ++++++ .../GeomBuilder/src/MagGeoBuilderFromDDD.cc | 13 +- .../GeomBuilder/src/MagGeoBuilderFromDDD.h | 51 +- MagneticField/GeomBuilder/src/bLayer.cc | 38 +- MagneticField/GeomBuilder/src/bLayer.h | 63 +- MagneticField/GeomBuilder/src/bRod.cc | 32 +- MagneticField/GeomBuilder/src/bRod.h | 36 +- MagneticField/GeomBuilder/src/bSector.cc | 61 +- MagneticField/GeomBuilder/src/bSector.h | 47 +- MagneticField/GeomBuilder/src/bSlab.cc | 17 +- MagneticField/GeomBuilder/src/bSlab.h | 60 +- MagneticField/GeomBuilder/src/buildBox.icc | 58 +- MagneticField/GeomBuilder/src/buildCons.icc | 52 +- .../GeomBuilder/src/buildPseudoTrap.icc | 102 ++- MagneticField/GeomBuilder/src/buildTrap.icc | 105 ++- .../GeomBuilder/src/buildTruncTubs.icc | 50 +- MagneticField/GeomBuilder/src/buildTubs.icc | 30 +- MagneticField/GeomBuilder/src/eLayer.cc | 8 +- MagneticField/GeomBuilder/src/eLayer.h | 33 +- MagneticField/GeomBuilder/src/eSector.cc | 22 +- MagneticField/GeomBuilder/src/eSector.h | 36 +- .../GeomBuilder/src/printUniqueNames.cc | 29 + .../GeomBuilder/src/printUniqueNames.h | 12 + MagneticField/GeomBuilder/src/volumeHandle.cc | 250 +------ MagneticField/GeomBuilder/src/volumeHandle.h | 197 +----- .../GeomBuilder/test/{ => python}/mfwriter.py | 0 .../test/{ => python}/mfxmlwriter.py | 0 .../test/python/testMagGeometry.py | 93 +++ .../test/python/testMagGeometryOldDD.py | 102 +++ .../test/stubs/MagGeometryAnalyzer.cc | 3 +- .../test/stubs/MagGeometryExerciser.cc | 2 +- .../interface/IdealMagneticFieldRecord.h | 8 +- .../VolumeBasedEngine/interface/MagGeometry.h | 1 - .../VolumeBasedEngine/src/MagGeometry.cc | 44 +- 44 files changed, 2393 insertions(+), 891 deletions(-) create mode 100644 DetectorDescription/DDCMS/plugins/DDCompactViewMFESProducer.cc create mode 100644 MagneticField/GeomBuilder/plugins/dd4hep/MagGeoBuilder.cc create mode 100644 MagneticField/GeomBuilder/plugins/dd4hep/MagGeoBuilder.h create mode 100644 MagneticField/GeomBuilder/plugins/dd4hep/VolBasedMagFieldESProducerNewDD.cc create mode 100644 MagneticField/GeomBuilder/plugins/dd4hep/volumeHandle.cc create mode 100644 MagneticField/GeomBuilder/plugins/dd4hep/volumeHandle.h create mode 100644 MagneticField/GeomBuilder/src/BaseVolumeHandle.cc create mode 100644 MagneticField/GeomBuilder/src/BaseVolumeHandle.h create mode 100644 MagneticField/GeomBuilder/src/printUniqueNames.cc create mode 100644 MagneticField/GeomBuilder/src/printUniqueNames.h rename MagneticField/GeomBuilder/test/{ => python}/mfwriter.py (100%) rename MagneticField/GeomBuilder/test/{ => python}/mfxmlwriter.py (100%) create mode 100644 MagneticField/GeomBuilder/test/python/testMagGeometry.py create mode 100644 MagneticField/GeomBuilder/test/python/testMagGeometryOldDD.py diff --git a/DetectorDescription/DDCMS/data/cms-mf-geometry.xml b/DetectorDescription/DDCMS/data/cms-mf-geometry.xml index b692a525c41bc..5c122a9aa1ca4 100644 --- a/DetectorDescription/DDCMS/data/cms-mf-geometry.xml +++ b/DetectorDescription/DDCMS/data/cms-mf-geometry.xml @@ -26,6 +26,7 @@ + diff --git a/DetectorDescription/DDCMS/plugins/BuildFile.xml b/DetectorDescription/DDCMS/plugins/BuildFile.xml index 7157051605199..ee47d29dcb644 100644 --- a/DetectorDescription/DDCMS/plugins/BuildFile.xml +++ b/DetectorDescription/DDCMS/plugins/BuildFile.xml @@ -17,6 +17,7 @@ + diff --git a/DetectorDescription/DDCMS/plugins/DDCompactViewMFESProducer.cc b/DetectorDescription/DDCMS/plugins/DDCompactViewMFESProducer.cc new file mode 100644 index 0000000000000..8a1fda1825547 --- /dev/null +++ b/DetectorDescription/DDCMS/plugins/DDCompactViewMFESProducer.cc @@ -0,0 +1,70 @@ +// -*- C++ -*- +// +// Package: DetectorDescription/Core +// Class: DDCompactViewMFESProducer +// +/**\class DDCompactViewMFESProducer + + Description: Produce DDCompactView + + Implementation: + Allow users view a DDDetector as a legacy compact view +*/ +// +// Original Author: Ianna Osborne +// +// + +#include + +#include "FWCore/Framework/interface/ModuleFactory.h" +#include "FWCore/Framework/interface/ESProducer.h" +#include "FWCore/Framework/interface/ESHandle.h" + +#include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h" +#include "FWCore/ParameterSet/interface/ParameterSetDescription.h" +#include "DetectorDescription/DDCMS/interface/DDCompactView.h" +#include "Geometry/Records/interface/GeometryFileRcd.h" +#include "DetectorDescription/DDCMS/interface/DDDetector.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "DD4hep/Detector.h" + +using namespace std; +using namespace cms; + +class DDCompactViewMFESProducer : public edm::ESProducer { +public: + DDCompactViewMFESProducer(const edm::ParameterSet&); + ~DDCompactViewMFESProducer() override; + + using ReturnType = unique_ptr; + + static void fillDescriptions(edm::ConfigurationDescriptions&); + + ReturnType produce(const IdealMagneticFieldRecord&); + +private: + const string m_label; +}; + +DDCompactViewMFESProducer::DDCompactViewMFESProducer(const edm::ParameterSet& iConfig) + : m_label(iConfig.getParameter("appendToDataLabel")) { + setWhatProduced(this); +} + +DDCompactViewMFESProducer::~DDCompactViewMFESProducer() {} + +void DDCompactViewMFESProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + edm::ParameterSetDescription desc; + descriptions.addDefault(desc); +} + +DDCompactViewMFESProducer::ReturnType DDCompactViewMFESProducer::produce(const IdealMagneticFieldRecord& iRecord) { + edm::ESHandle det; + iRecord.getRecord().get(m_label, det); + + auto product = std::make_unique(*det); + return product; +} + +DEFINE_FWK_EVENTSETUP_MODULE(DDCompactViewMFESProducer); diff --git a/MagneticField/GeomBuilder/plugins/BuildFile.xml b/MagneticField/GeomBuilder/plugins/BuildFile.xml index fa5245600e2c0..aa7c983ed533c 100644 --- a/MagneticField/GeomBuilder/plugins/BuildFile.xml +++ b/MagneticField/GeomBuilder/plugins/BuildFile.xml @@ -7,3 +7,10 @@ + + + + + + + diff --git a/MagneticField/GeomBuilder/plugins/dd4hep/MagGeoBuilder.cc b/MagneticField/GeomBuilder/plugins/dd4hep/MagGeoBuilder.cc new file mode 100644 index 0000000000000..a73d1e5fc9288 --- /dev/null +++ b/MagneticField/GeomBuilder/plugins/dd4hep/MagGeoBuilder.cc @@ -0,0 +1,631 @@ +/* + * See header file for a description of this class. + * + */ + +#include "MagneticField/GeomBuilder/plugins/dd4hep/MagGeoBuilder.h" +#include "MagneticField/GeomBuilder/src/bLayer.h" +#include "MagneticField/GeomBuilder/src/eSector.h" +#include "MagneticField/GeomBuilder/src/FakeInterpolator.h" + +#include "MagneticField/Layers/interface/MagBLayer.h" +#include "MagneticField/Layers/interface/MagESector.h" + +#include "FWCore/ParameterSet/interface/FileInPath.h" + +#include "DetectorDescription/DDCMS/interface/DDFilteredView.h" + +#include "Utilities/BinningTools/interface/ClusterizingHistogram.h" + +#include "MagneticField/Interpolation/interface/MagProviderInterpol.h" +#include "MagneticField/Interpolation/interface/MFGridFactory.h" +#include "MagneticField/Interpolation/interface/MFGrid.h" + +#include "MagneticField/VolumeGeometry/interface/MagVolume6Faces.h" +#include "MagneticField/VolumeGeometry/interface/MagExceptions.h" +#include "MagneticField/Layers/interface/MagVerbosity.h" + +#include "DataFormats/Math/interface/deltaPhi.h" + +#include "FWCore/Utilities/interface/Exception.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "Utilities/General/interface/precomputed_value_sort.h" +#include "DD4hep/Detector.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace std; +using namespace magneticfield; +using namespace cms; +using namespace edm; +using namespace dd4hep; +using namespace angle_units::operators; + +MagGeoBuilder::MagGeoBuilder(string tableSet, int geometryVersion, bool debug) + : tableSet_(tableSet), geometryVersion_(geometryVersion), theGridFiles_(nullptr), debug_(debug) { + if (debug_) + LogTrace("MagGeoBuilder") << "Constructing a MagGeoBuilder" << endl; +} + +MagGeoBuilder::~MagGeoBuilder() { + for (handles::const_iterator i = bVolumes_.begin(); i != bVolumes_.end(); ++i) { + delete (*i); + } + + for (handles::const_iterator i = eVolumes_.begin(); i != eVolumes_.end(); ++i) { + delete (*i); + } +} + +void MagGeoBuilder::summary(handles& volumes) { + // The final countdown. + int ivolumes = volumes.size(); // number of volumes + int isurfaces = ivolumes * 6; // number of individual surfaces + int iassigned = 0; // How many have been assigned + int iunique = 0; // number of unique surfaces + int iref_ass = 0; + int iref_nass = 0; + + set ptrs; + + handles::const_iterator first = volumes.begin(); + handles::const_iterator last = volumes.end(); + + for (handles::const_iterator i = first; i != last; ++i) { + DDSolidShape theShape = (*i)->shape(); + if (theShape == DDSolidShape::ddbox || theShape == DDSolidShape::ddcons || theShape == DDSolidShape::ddtrap || + theShape == DDSolidShape::ddtubs) { + for (int side = 0; side < 6; ++side) { + int references = (*i)->references(side); + if ((*i)->isPlaneMatched(side)) { + ++iassigned; + bool firstOcc = (ptrs.insert(&((*i)->surface(side)))).second; + if (firstOcc) + iref_ass += references; + if (references < 2) { + LogTrace("MagGeoBuilder") << "*** Only 1 ref, vol: " << (*i)->volumeno << " # " << (*i)->copyno + << " side: " << side << endl; + } + } else { + iref_nass += references; + if (references > 1) { + LogTrace("MagGeoBuilder") << "*** Ref_nass >1 " << endl; + } + } + } + } // end if theShape + } // end for + iunique = ptrs.size(); + + LogTrace("MagGeoBuilder") << " volumes " << ivolumes << endl + << " surfaces " << isurfaces << endl + << " assigned " << iassigned << endl + << " unique " << iunique << endl + << " iref_ass " << iref_ass << endl + << " iref_nass " << iref_nass << endl; +} + +void MagGeoBuilder::build(const DDDetector* det) { + Volume top = det->worldVolume(); + DDFilteredView fv(det, top); + if (fv.next(0) == false) { + LogError("MagGeoBuilder") << "Filtered view is empty. Cannot build."; + return; + } + + // The actual field interpolators + map bInterpolators; + map eInterpolators; + + // Counter of different volumes + int bVolCount = 0; + int eVolCount = 0; + + const string magfStr{"MAGF"}; + if (fv.name() != magfStr) { + std::string topNodeName(fv.name()); + LogTrace("MagGeoBuilder") << "Filtered view top node name is " << topNodeName << "."; + + //see if one of the children is MAGF + bool doSubDets = fv.next(0); + + bool go = true; + while (go && doSubDets) { + LogTrace("MagGeoBuilder") << "Next node name is " << fv.name() << "."; + if (fv.name() == magfStr) + break; + else + go = fv.next(0); + } + if (!go) { + throw cms::Exception("NoMAGF") + << " Neither the top node, nor any child node of the filtered view is \"MAGF\" but the top node is instead \"" + << topNodeName << "\""; + } + } + + // Loop over MAGF volumes and create volumeHandles. + bool doSubDets = fv.next(0); + if (doSubDets == false) { + LogError("MagGeoBuilder") << "Filtered view has no node. Cannot build."; + return; + } + while (doSubDets) { + string name = fv.volume().volume().name(); + LogTrace("MagGeoBuilder") << "Name: " << name; + + bool expand = false; + volumeHandle* v = new volumeHandle(fv, expand, debug_); + + if (theGridFiles_ != nullptr) { + int key = (v->volumeno) * 100 + v->copyno; + TableFileMap::const_iterator itable = theGridFiles_->find(key); + if (itable == theGridFiles_->end()) { + key = (v->volumeno) * 100; + itable = theGridFiles_->find(key); + } + + if (itable != theGridFiles_->end()) { + string magFile = (*itable).second.first; + stringstream conv; + string svol, ssec; + conv << setfill('0') << setw(3) << v->volumeno << " " << setw(2) + << v->copyno; // volume assumed to have 0s padding to 3 digits; sector assumed to have 0s padding to 2 digits + conv >> svol >> ssec; + boost::replace_all(magFile, "[v]", svol); + boost::replace_all(magFile, "[s]", ssec); + int masterSector = (*itable).second.second; + if (masterSector == 0) + masterSector = v->copyno; + v->magFile = magFile; + v->masterSector = masterSector; + } else { + edm::LogError("MagGeoBuilderbuild") << "ERROR: no table spec found for V " << v->volumeno << ":" << v->copyno; + } + } + + // Select volumes, build volume handles. + float Z = v->center().z(); + float R = v->center().perp(); + LogTrace("MagGeoBuilder") << " Vol R and Z values determine barrel or endcap. R = " << R << ", Z = " << Z; + + // v 85l: Barrel is everything up to |Z| = 661.0, excluding + // volume #7, centered at 6477.5 + // v 1103l: same numbers work fine. #16 instead of #7, same coords; + // see comment below for V6,7 + //ASSUMPTION: no misalignment is applied to mag volumes. + //FIXME: implement barrel/endcap flags as DDD SpecPars. + if ((fabs(Z) < 647. || (R > 350. && fabs(Z) < 662.)) && + !(fabs(Z) > 480 && R < 172) // in 1103l we place V_6 and V_7 in the + // endcaps to preserve nice layer structure + // in the barrel. This does not hurt in v85l + // where there is a single V1 + ) { // Barrel + if (debug_) + LogTrace("MagGeoBuilder") << " (Barrel)" << endl; + bVolumes_.push_back(v); + + // Build the interpolator of the "master" volume (the one which is + // not replicated in phi) + // ASSUMPTION: copyno == sector. + if (v->copyno == v->masterSector) { + buildInterpolator(v, bInterpolators); + ++bVolCount; + } + } else { // Endcaps + if (debug_) + LogTrace("MagGeoBuilder") << " (Endcaps)" << endl; + eVolumes_.push_back(v); + if (v->copyno == v->masterSector) { + buildInterpolator(v, eInterpolators); + ++eVolCount; + } + } + doSubDets = fv.next(0); // end of loop over MAGF + } + + if (debug_) { + LogTrace("MagGeoBuilder") << "Number of volumes (barrel): " << bVolumes_.size() << endl + << " Number of volumes (endcap): " << eVolumes_.size() << endl; + LogTrace("MagGeoBuilder") << "**********************************************************" << endl; + } + + // Now all volumeHandles are there, and parameters for each of the planes + // are calculated. + + //---------------------------------------------------------------------- + // Print summary information + + if (debug_) { + LogTrace("MagGeoBuilder") << "-----------------------" << endl; + LogTrace("MagGeoBuilder") << "SUMMARY: Barrel " << endl; + summary(bVolumes_); + + LogTrace("MagGeoBuilder") << endl << "SUMMARY: Endcaps " << endl; + summary(eVolumes_); + LogTrace("MagGeoBuilder") << "-----------------------" << endl; + } + + //---------------------------------------------------------------------- + // Find barrel layers. + + if (bVolumes_.empty()) { + LogError("MagGeoBuilder") << "Error: Barrel volumes are missing. Terminating build."; + return; + } + vector layers; // the barrel layers + precomputed_value_sort(bVolumes_.begin(), bVolumes_.end(), ExtractRN()); + + // Find the layers (in R) + const float resolution = 1.; // cm + float rmin = bVolumes_.front()->RN() - resolution; + float rmax = bVolumes_.back()->RN() + resolution; + ClusterizingHistogram hisR(int((rmax - rmin) / resolution) + 1, rmin, rmax); + + if (debug_) + LogTrace("MagGeoBuilder") << " R layers: " << rmin << " " << rmax << endl; + + handles::const_iterator first = bVolumes_.begin(); + handles::const_iterator last = bVolumes_.end(); + + for (handles::const_iterator i = first; i != last; ++i) { + hisR.fill((*i)->RN()); + } + vector rClust = hisR.clusterize(resolution); + + handles::const_iterator ringStart = first; + handles::const_iterator separ = first; + + for (unsigned int i = 0; i < rClust.size() - 1; ++i) { + if (debug_) + LogTrace("MagGeoBuilder") << " Layer at RN = " << rClust[i]; + float rSepar = (rClust[i] + rClust[i + 1]) / 2.f; + while ((*separ)->RN() < rSepar) + ++separ; + + bLayer thislayer(ringStart, separ, debug_); + layers.push_back(thislayer); + ringStart = separ; + } + { + if (debug_) + LogTrace("MagGeoBuilder") << " Layer at RN = " << rClust.back(); + bLayer thislayer(separ, last, debug_); + layers.push_back(thislayer); + } + + if (debug_) + LogTrace("MagGeoBuilder") << "Barrel: Found " << rClust.size() << " clusters in R, " << layers.size() + << " layers " << endl + << endl; + + //---------------------------------------------------------------------- + // Find endcap sectors + vector sectors; // the endcap sectors + + // Find the number of sectors (should be 12 or 24 depending on the geometry model) + constexpr float phireso = 0.05; // rad + constexpr int twoPiOverPhiReso = static_cast(2._pi / phireso) + 1; + ClusterizingHistogram hisPhi(twoPiOverPhiReso, -1._pi, 1._pi); + + for (handles::const_iterator i = eVolumes_.begin(); i != eVolumes_.end(); ++i) { + hisPhi.fill((*i)->minPhi()); + } + vector phiClust = hisPhi.clusterize(phireso); + int nESectors = phiClust.size(); + if (nESectors <= 0) { + LogError("MagGeoBuilder") << "ERROR: Endcap sectors are missing. Terminating build."; + return; + } + if (debug_ && (nESectors % 12) != 0) { + LogTrace("MagGeoBuilder") << "ERROR: unexpected # of endcap sectors: " << nESectors; + } + + //Sort in phi + precomputed_value_sort(eVolumes_.begin(), eVolumes_.end(), ExtractPhi()); + + // Handle the -pi/pi boundary: volumes crossing it could be half at the begin and half at end of the sorted list. + // So, check if any of the volumes that should belong to the first bin (at -phi) are at the end of the list: + float lastBinPhi = phiClust.back(); + handles::reverse_iterator ri = eVolumes_.rbegin(); + while ((*ri)->center().phi() > lastBinPhi) { + ++ri; + } + if (ri != eVolumes_.rbegin()) { + // ri points to the first element that is within the last bin. + // We need to move the following element (ie ri.base()) to the beginning of the list, + handles::iterator newbeg = ri.base(); + rotate(eVolumes_.begin(), newbeg, eVolumes_.end()); + } + + //Group volumes in sectors + int offset = eVolumes_.size() / nESectors; + for (int i = 0; i < nESectors; ++i) { + if (debug_) { + LogTrace("MagGeoBuilder") << " Sector at phi = " << (*(eVolumes_.begin() + ((i)*offset)))->center().phi() + << endl; + // Additional x-check: sectors are expected to be made by volumes with the same copyno + int secCopyNo = -1; + for (handles::const_iterator iv = eVolumes_.begin() + ((i)*offset); iv != eVolumes_.begin() + ((i + 1) * offset); + ++iv) { + if (secCopyNo >= 0 && (*iv)->copyno != secCopyNo) + LogTrace("MagGeoBuilder") << "ERROR: volume copyno " << (*iv)->name << ":" << (*iv)->copyno + << " differs from others in same sectors with copyno = " << secCopyNo << endl; + secCopyNo = (*iv)->copyno; + } + } + + sectors.push_back(eSector(eVolumes_.begin() + ((i)*offset), eVolumes_.begin() + ((i + 1) * offset), debug_)); + } + + if (debug_) + LogTrace("MagGeoBuilder") << "Endcap: Found " << sectors.size() << " sectors " << endl; + + //---------------------------------------------------------------------- + // Build MagVolumes and the MagGeometry hierarchy. + + //--- Barrel + + // Build MagVolumes and associate interpolators to them + buildMagVolumes(bVolumes_, bInterpolators); + + // Build MagBLayers + for (vector::const_iterator ilay = layers.begin(); ilay != layers.end(); ++ilay) { + mBLayers_.push_back((*ilay).buildMagBLayer()); + } + + if (debug_) { + LogTrace("MagGeoBuilder") << "*** BARREL ********************************************" << endl + << "Number of different volumes = " << bVolCount << endl + << "Number of interpolators built = " << bInterpolators.size() << endl + << "Number of MagBLayers built = " << mBLayers_.size() << endl; + + testInside(bVolumes_); // FIXME: all volumes should be checked in one go. + } + + //--- Endcap + // Build MagVolumes and associate interpolators to them + buildMagVolumes(eVolumes_, eInterpolators); + + // Build the MagESectors + for (vector::const_iterator isec = sectors.begin(); isec != sectors.end(); ++isec) { + mESectors_.push_back((*isec).buildMagESector()); + } + + if (debug_) { + LogTrace("MagGeoBuilder") << "*** ENDCAP ********************************************" << endl + << "Number of different volumes = " << eVolCount << endl + << "Number of interpolators built = " << eInterpolators.size() << endl + << "Number of MagESector built = " << mESectors_.size() << endl; + + testInside(eVolumes_); // FIXME: all volumes should be checked in one go. + } +} + +void MagGeoBuilder::buildMagVolumes(const handles& volumes, map& interpolators) { + // Build all MagVolumes setting the MagProviderInterpol + for (handles::const_iterator vol = volumes.begin(); vol != volumes.end(); ++vol) { + const MagProviderInterpol* mp = nullptr; + if (interpolators.find((*vol)->magFile) != interpolators.end()) { + mp = interpolators[(*vol)->magFile]; + } else { + edm::LogError("MagGeoBuilder|buildMagVolumes") + << "No interpolator found for file " << (*vol)->magFile << " vol: " << (*vol)->volumeno << "\n" + << interpolators.size() << endl; + } + + // Search for [volume,sector] in the list of scaling factors; sector = 0 handled as wildcard + // ASSUMPTION: copyno == sector. + int key = ((*vol)->volumeno) * 100 + (*vol)->copyno; + map::const_iterator isf = theScalingFactors_.find(key); + if (isf == theScalingFactors_.end()) { + key = ((*vol)->volumeno) * 100; + isf = theScalingFactors_.find(key); + } + + double sf = 1.; + if (isf != theScalingFactors_.end()) { + sf = (*isf).second; + + LogTrace("MagGeoBuilder|buildMagVolumes") + << "Applying scaling factor " << sf << " to " << (*vol)->volumeno << "[" << (*vol)->copyno << "] (key:" << key + << ")" << endl; + } + + const GloballyPositioned* gpos = (*vol)->placement(); + (*vol)->magVolume = new MagVolume6Faces(gpos->position(), gpos->rotation(), (*vol)->sides(), mp, sf); + + if ((*vol)->copyno == (*vol)->masterSector) { + (*vol)->magVolume->ownsFieldProvider(true); + } + + (*vol)->magVolume->setIsIron((*vol)->isIron()); + + // The name and sector of the volume are saved for debug purposes only. They may be removed at some point... + (*vol)->magVolume->volumeNo = (*vol)->volumeno; + (*vol)->magVolume->copyno = (*vol)->copyno; + } +} + +void MagGeoBuilder::buildInterpolator(const volumeHandle* vol, map& interpolators) { + // Phi of the master sector + double masterSectorPhi = (vol->masterSector - 1) * 1._pi / 6.; + + if (debug_) { + LogTrace("MagGeoBuilder") << "Building interpolator from " << vol->volumeno << " copyno " << vol->copyno + << " at " << vol->center() + << " phi: " << static_cast(vol->center().phi()) / 1._pi + << " pi, file: " << vol->magFile << " master: " << vol->masterSector << endl; + + double delta = std::abs(vol->center().phi() - masterSectorPhi); + if (delta > (1._pi / 9.)) { + LogTrace("MagGeoBuilder") << "***WARNING wrong sector? Vol delta from master sector is " << delta / 1._pi + << " pi"; + } + } + + if (tableSet_ == "fake" || vol->magFile == "fake") { + interpolators[vol->magFile] = new magneticfield::FakeInterpolator(); + return; + } + + string fullPath; + + try { + edm::FileInPath mydata("MagneticField/Interpolation/data/" + tableSet_ + "/" + vol->magFile); + fullPath = mydata.fullPath(); + } catch (edm::Exception& exc) { + cerr << "MagGeoBuilder: exception in reading table; " << exc.what() << endl; + if (!debug_) + throw; + return; + } + + try { + if (vol->toExpand()) { + //FIXME: see discussion on mergeCylinders above. + // interpolators[vol->magFile] = + // MFGridFactory::build( fullPath, *(vol->placement()), vol->minPhi(), vol->maxPhi()); + } else { + // If the table is in "local" coordinates, must create a reference + // frame that is appropriately rotated along the CMS Z axis. + + GloballyPositioned rf = *(vol->placement()); + + if (vol->masterSector != 1) { + typedef Basic3DVector Vector; + + GloballyPositioned::RotationType rot(Vector(0, 0, 1), -masterSectorPhi); + Vector vpos(vol->placement()->position()); + + rf = GloballyPositioned(GloballyPositioned::PositionType(rot.multiplyInverse(vpos)), + vol->placement()->rotation() * rot); + } + + interpolators[vol->magFile] = MFGridFactory::build(fullPath, rf); + } + } catch (MagException& exc) { + LogTrace("MagGeoBuilder") << exc.what() << endl; + interpolators.erase(vol->magFile); + if (!debug_) + throw; + return; + } + + if (debug_) { + // Check that all grid points of the interpolator are inside the volume. + const MagVolume6Faces tempVolume( + vol->placement()->position(), vol->placement()->rotation(), vol->sides(), interpolators[vol->magFile]); + + const MFGrid* grid = dynamic_cast(interpolators[vol->magFile]); + if (grid != nullptr) { + Dimensions sizes = grid->dimensions(); + LogTrace("MagGeoBuilder") << "Grid has 3 dimensions " + << " number of nodes is " << sizes.w << " " << sizes.h << " " << sizes.d << endl; + + const double tolerance = 0.03; + + size_t dumpCount = 0; + for (int j = 0; j < sizes.h; j++) { + for (int k = 0; k < sizes.d; k++) { + for (int i = 0; i < sizes.w; i++) { + MFGrid::LocalPoint lp = grid->nodePosition(i, j, k); + if (!tempVolume.inside(lp, tolerance)) { + if (++dumpCount < 2) { + MFGrid::GlobalPoint gp = tempVolume.toGlobal(lp); + LogTrace("MagGeoBuilder") << "GRID ERROR: " << i << " " << j << " " << k << " local: " << lp + << " global: " << gp << " R= " << gp.perp() << " phi=" << gp.phi() << endl; + } + } + } + } + } + + LogTrace("MagGeoBuilder") << "Volume:" << vol->volumeno + << " : Number of grid points outside the MagVolume: " << dumpCount << "/" + << sizes.w * sizes.h * sizes.d << endl; + } + } +} + +void MagGeoBuilder::testInside(handles& volumes) { + // test inside() for all volumes. + LogTrace("MagGeoBuilder") << "--------------------------------------------------" << endl; + LogTrace("MagGeoBuilder") << " inside(center) test" << endl; + for (handles::const_iterator vol = volumes.begin(); vol != volumes.end(); ++vol) { + for (handles::const_iterator i = volumes.begin(); i != volumes.end(); ++i) { + if ((*i) == (*vol)) + continue; + //if ((*i)->magVolume == 0) continue; + if ((*i)->magVolume->inside((*vol)->center())) { + LogTrace("MagGeoBuilder") << "*** ERROR: center of V " << (*vol)->volumeno << ":" << (*vol)->copyno + << " is inside V " << (*i)->volumeno << ":" << (*i)->copyno << endl; + } + } + + if ((*vol)->magVolume->inside((*vol)->center())) { + LogTrace("MagGeoBuilder") << "V " << (*vol)->volumeno << ":" << (*vol)->copyno << " OK " << endl; + } else { + LogTrace("MagGeoBuilder") << "*** ERROR: center of volume is not inside it, " << (*vol)->volumeno << ":" + << (*vol)->copyno; + } + } + LogTrace("MagGeoBuilder") << "--------------------------------------------------" << endl; +} + +vector MagGeoBuilder::barrelLayers() const { return mBLayers_; } + +vector MagGeoBuilder::endcapSectors() const { return mESectors_; } + +vector MagGeoBuilder::barrelVolumes() const { + vector v; + v.reserve(bVolumes_.size()); + for (handles::const_iterator i = bVolumes_.begin(); i != bVolumes_.end(); ++i) { + v.push_back((*i)->magVolume); + } + return v; +} + +vector MagGeoBuilder::endcapVolumes() const { + vector v; + v.reserve(eVolumes_.size()); + for (handles::const_iterator i = eVolumes_.begin(); i != eVolumes_.end(); ++i) { + v.push_back((*i)->magVolume); + } + return v; +} + +float MagGeoBuilder::maxR() const { + //FIXME: should get it from the actual geometry + return 900.; +} + +float MagGeoBuilder::maxZ() const { + //FIXME: should get it from the actual geometry + if (geometryVersion_ >= 160812) + return 2400.; + else if (geometryVersion_ >= 120812) + return 2000.; + else + return 1600.; +} + +void MagGeoBuilder::setScaling(const std::vector& keys, const std::vector& values) { + if (keys.size() != values.size()) { + throw cms::Exception("InvalidParameter") + << "Invalid field scaling parameters 'scalingVolumes' and 'scalingFactors' "; + } + for (unsigned int i = 0; i < keys.size(); ++i) { + theScalingFactors_[keys[i]] = values[i]; + } +} + +void MagGeoBuilder::setGridFiles(const TableFileMap& gridFiles) { theGridFiles_ = &gridFiles; } diff --git a/MagneticField/GeomBuilder/plugins/dd4hep/MagGeoBuilder.h b/MagneticField/GeomBuilder/plugins/dd4hep/MagGeoBuilder.h new file mode 100644 index 0000000000000..af1dec3bef5e8 --- /dev/null +++ b/MagneticField/GeomBuilder/plugins/dd4hep/MagGeoBuilder.h @@ -0,0 +1,98 @@ +#ifndef MagneticField_GeomBuilder_MagGeoBuilder_H +#define MagneticField_GeomBuilder_MagGeoBuilder_H + +/** \class MagGeoBuilder + * Parse the XML magnetic geometry, build individual volumes and match their + * shared surfaces. Build MagVolume6Faces and organise them in a hierarchical + * structure. Build MagGeometry out of it. + * + * \author N. Amapane - INFN Torino + */ + +#include "CondFormats/MFObjects/interface/MagFieldConfig.h" +#include "DataFormats/GeometrySurface/interface/ReferenceCounted.h" +#include "DetectorDescription/DDCMS/interface/DDSpecParRegistry.h" +#include "DetectorDescription/DDCMS/interface/DDDetector.h" +#include "Geometry/MuonNumbering/interface/DD4hep_DTNumberingScheme.h" +#include "MagneticField/Interpolation/interface/MagProviderInterpol.h" +#include "MagneticField/GeomBuilder/plugins/dd4hep/volumeHandle.h" + +#include +#include +#include +#include + +namespace dd4hep { + class Detector; +} + +class Surface; +class MagBLayer; +class MagESector; +class MagVolume6Faces; + +namespace cms { + + class MagGeoBuilder { + public: + using handles = magneticfield::handles; + + MagGeoBuilder(std::string tableSet, int geometryVersion, bool debug = false); + + ~MagGeoBuilder(); + + // Build the geometry. + void build(const DDDetector* det); + + /// Set scaling factors for individual volumes. + /// "keys" is a vector of 100*volume number + sector (sector 0 = all sectors) + /// "values" are the corresponding scaling factors + void setScaling(const std::vector& keys, const std::vector& values); + + void setGridFiles(const magneticfield::TableFileMap& gridFiles); + + /// Get barrel layers + std::vector barrelLayers() const; + + /// Get endcap layers + std::vector endcapSectors() const; + + float maxR() const; + + float maxZ() const; + + std::vector barrelVolumes() const; + std::vector endcapVolumes() const; + + private: + typedef ConstReferenceCountingPointer RCPS; + + // Build interpolator for the volume with "correct" rotation + void buildInterpolator(const magneticfield::volumeHandle* vol, + std::map& interpolators); + + // Build all MagVolumes setting the MagProviderInterpol + void buildMagVolumes(const handles& volumes, std::map& interpolators); + + // Print checksums for surfaces. + void summary(handles& volumes); + + // Perform simple sanity checks + void testInside(handles& volumes); + + handles bVolumes_; // the barrel volumes. + handles eVolumes_; // the endcap volumes. + + std::vector mBLayers_; // Finally built barrel geometry + std::vector mESectors_; // Finally built barrel geometry + + std::string tableSet_; // Version of the data files to be used + int geometryVersion_; // Version of MF geometry + + std::map theScalingFactors_; + const magneticfield::TableFileMap* theGridFiles_; // Non-owned pointer assumed to be valid until build() is called + + const bool debug_; + }; +} // namespace cms +#endif diff --git a/MagneticField/GeomBuilder/plugins/dd4hep/VolBasedMagFieldESProducerNewDD.cc b/MagneticField/GeomBuilder/plugins/dd4hep/VolBasedMagFieldESProducerNewDD.cc new file mode 100644 index 0000000000000..2af5ba1bce488 --- /dev/null +++ b/MagneticField/GeomBuilder/plugins/dd4hep/VolBasedMagFieldESProducerNewDD.cc @@ -0,0 +1,114 @@ +/** \class VolBasedMagFieldESProducerNewDD + * + * Producer for the VolumeBasedMagneticField. + * + */ + +#include "FWCore/Framework/interface/ESTransientHandle.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/Framework/interface/EventSetupRecordIntervalFinder.h" +#include "FWCore/Framework/interface/ESProducer.h" +#include "FWCore/Framework/interface/ModuleFactory.h" + +#include "MagneticField/VolumeBasedEngine/interface/VolumeBasedMagneticField.h" +#include "MagneticField/Records/interface/IdealMagneticFieldRecord.h" +#include "MagneticField/Engine/interface/MagneticField.h" + +#include "MagneticField/GeomBuilder/plugins/dd4hep/MagGeoBuilder.h" +#include "CondFormats/MFObjects/interface/MagFieldConfig.h" +#include "Geometry/Records/interface/GeometryFileRcd.h" +#include "DetectorDescription/DDCMS/interface/BenchmarkGrd.h" +#include "DetectorDescription/DDCMS/interface/DDCompactView.h" +#include "DetectorDescription/DDCMS/interface/DDDetector.h" + +#include +#include +#include + +using namespace cms; +using namespace std; +using namespace magneticfield; + +namespace magneticfield { + class VolBasedMagFieldESProducerNewDD : public edm::ESProducer { + public: + VolBasedMagFieldESProducerNewDD(const edm::ParameterSet& iConfig); + + // forbid copy ctor and assignment op. + VolBasedMagFieldESProducerNewDD(const VolBasedMagFieldESProducerNewDD&) = delete; + const VolBasedMagFieldESProducerNewDD& operator=(const VolBasedMagFieldESProducerNewDD&) = delete; + + std::unique_ptr produce(const IdealMagneticFieldRecord& iRecord); + + private: + edm::ParameterSet pset_; + const bool debug_; + const bool useParametrizedTrackerField_; + const MagFieldConfig conf_; + const std::string version_; + edm::ESGetToken paramFieldToken_; + edm::ESGetToken cpvToken_; + const edm::ESInputTag tag_; + }; +} // namespace magneticfield + +VolBasedMagFieldESProducerNewDD::VolBasedMagFieldESProducerNewDD(const edm::ParameterSet& iConfig) + : pset_{iConfig}, + debug_{iConfig.getUntrackedParameter("debugBuilder", false)}, + useParametrizedTrackerField_{iConfig.getParameter("useParametrizedTrackerField")}, + conf_{iConfig, debug_}, + version_{iConfig.getParameter("version")}, + tag_{iConfig.getParameter("DDDetector")} { + // LogInfo used because LogDebug messages don't appear even when fully enabled. + edm::LogInfo("VolBasedMagFieldESProducerNewDD") << "info:Constructing a VolBasedMagFieldESProducerNewDD" << endl; + + auto cc = setWhatProduced(this, iConfig.getUntrackedParameter("label", "")); + cc.setConsumes(cpvToken_, edm::ESInputTag{"", "magfield"}); + if (useParametrizedTrackerField_) { + cc.setConsumes(paramFieldToken_, edm::ESInputTag{"", iConfig.getParameter("paramLabel")}); + } +} + +// ------------ method called to produce the data ------------ +std::unique_ptr VolBasedMagFieldESProducerNewDD::produce(const IdealMagneticFieldRecord& iRecord) { + if (debug_) { + edm::LogInfo("VolBasedMagFieldESProducerNewDD") << "VolBasedMagFieldESProducerNewDD::produce() " << version_; + } + + MagGeoBuilder builder(conf_.version, conf_.geometryVersion, debug_); + + // Set scaling factors + if (!conf_.keys.empty()) { + builder.setScaling(conf_.keys, conf_.values); + } + + // Set specification for the grid tables to be used. + if (!conf_.gridFiles.empty()) { + builder.setGridFiles(conf_.gridFiles); + } + + auto cpv = iRecord.getTransientHandle(cpvToken_); + const DDCompactView* cpvPtr = cpv.product(); + const DDDetector* det = cpvPtr->detector(); + builder.build(det); + edm::LogInfo("VolBasedMagFieldESProducerNewDD") << "produce() finished build"; + + // Get slave field (from ES) + const MagneticField* paramField = nullptr; + if (useParametrizedTrackerField_) { + edm::LogInfo("VolBasedMagFieldESProducerNewDD") << "Getting MF for parametrized field"; + paramField = &iRecord.get(paramFieldToken_); + } + return std::make_unique(conf_.geometryVersion, + builder.barrelLayers(), + builder.endcapSectors(), + builder.barrelVolumes(), + builder.endcapVolumes(), + builder.maxR(), + builder.maxZ(), + paramField, + false); +} + +DEFINE_FWK_EVENTSETUP_MODULE(VolBasedMagFieldESProducerNewDD); diff --git a/MagneticField/GeomBuilder/plugins/dd4hep/volumeHandle.cc b/MagneticField/GeomBuilder/plugins/dd4hep/volumeHandle.cc new file mode 100644 index 0000000000000..2846b54d421eb --- /dev/null +++ b/MagneticField/GeomBuilder/plugins/dd4hep/volumeHandle.cc @@ -0,0 +1,207 @@ +// #include "Utilities/Configuration/interface/Architecture.h" + +/* + * See header file for a description of this class. + * + * \author N. Amapane - INFN Torino + */ + +#include "MagneticField/GeomBuilder/plugins/dd4hep/volumeHandle.h" + +#include "DataFormats/GeometrySurface/interface/Plane.h" +#include "DataFormats/GeometrySurface/interface/Cylinder.h" +#include "DataFormats/GeometrySurface/interface/Cone.h" +#include "DataFormats/GeometryVector/interface/CoordinateSets.h" +#include "DataFormats/Math/interface/GeantUnits.h" +#include "DataFormats/Math/interface/Vector3D.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" + +#include "DetectorDescription/DDCMS/interface/DDShapes.h" + +#include "MagneticField/Layers/interface/MagVerbosity.h" + +#include +#include +#include +#include +#include + +using namespace SurfaceOrientation; +using namespace std; +using namespace magneticfield; +using namespace cms; +using namespace cms::dd; +using namespace edm; + +volumeHandle::volumeHandle(const DDFilteredView &fv, bool expand2Pi, bool debugVal) + : BaseVolumeHandle(debugVal), shape_(getCurrentShape(fv)), solid(fv) { + name = fv.name(); + copyno = fv.copyNum(); + expand = expand2Pi; + const Double_t *const transArray = fv.trans(); + center_ = GlobalPoint(transArray[0], transArray[1], transArray[2]); + + // ASSUMPTION: volume names ends with "_NUM" where NUM is the volume number + string volName = name; + volName.erase(0, volName.rfind('_') + 1); + volumeno = boost::lexical_cast(volName); + + for (int i = 0; i < 6; ++i) { + isAssigned[i] = false; + } + referencePlane(fv); + switch (shape_) { + case DDSolidShape::ddbox: + buildBox(); + break; + case DDSolidShape::ddtrap: + buildTrap(); + break; + case DDSolidShape::ddcons: + buildCons(); + break; + case DDSolidShape::ddtubs: + buildTubs(); + break; + case DDSolidShape::ddtrunctubs: + buildTruncTubs(); + break; + default: + LogError("magneticfield::volumeHandle") + << "ctor: Unexpected shape_ # " << static_cast(shape_) << " for vol " << name; + } + + // Get material for this volume + if (fv.materialName() == "Iron") + isIronFlag = true; + + if (debug) { + LogDebug("magneticfield::volumeHandle") << " RMin = " << theRMin; + LogDebug("magneticfield::volumeHandle") << " RMax = " << theRMax; + + if (theRMin < 0 || theRN < theRMin || theRMax < theRN) + LogDebug("magneticfield::volumeHandle") << "*** WARNING: wrong RMin/RN/RMax"; + + LogDebug("magneticfield::volumeHandle") + << "Summary: " << name << " " << copyno << " shape = " << shape_ << " trasl " << center() << " R " + << center().perp() << " phi " << center().phi() << " magFile " << magFile << " Material= " << fv.materialName() + << " isIron= " << isIronFlag << " masterSector= " << masterSector << std::endl; + + LogDebug("magneticfield::volumeHandle") << " Orientation of surfaces:"; + std::string sideName[3] = {"positiveSide", "negativeSide", "onSurface"}; + for (int i = 0; i < 6; ++i) { + if (surfaces[i] != nullptr) + LogDebug("magneticfield::volumeHandle") << " " << i << ":" << sideName[surfaces[i]->side(center_, 0.3)]; + } + } +} + +void volumeHandle::referencePlane(const DDFilteredView &fv) { + // The refPlane is the "main plane" for the solid. It corresponds to the + // x,y plane in the DDD local frame, and defines a frame where the local + // coordinates are the same as in DDD. + // In the geometry version 85l_030919, this plane is normal to the + // beam line for all volumes but pseudotraps, so that global R is along Y, + // global phi is along -X and global Z along Z: + // + // Global(for vol at pi/2) Local + // +R (+Y) +Y + // +phi(-X) -X + // +Z +Z + // + // For pseudotraps the refPlane is parallel to beam line and global R is + // along Z, global phi is along +-X and and global Z along Y: + // + // Global(for vol at pi/2) Local + // +R (+Y) +Z + // +phi(-X) +X + // +Z +Y + // + // Note that the frame is centered in the DDD volume center, which is + // inside the volume for DDD boxes and (pesudo)trapezoids, on the beam line + // for tubs, cons and trunctubs. + + // In geometry version 1103l, trapezoids have X and Z in the opposite direction + // than the above. Boxes are either oriented as described above or in some case + // have opposite direction for Y and X. + + // The global position + Surface::PositionType &posResult = center_; + + // The reference plane rotation + math::XYZVector x, y, z; + dd4hep::Rotation3D refRot; + fv.rot(refRot); + refRot.GetComponents(x, y, z); + if (debug) { + if (x.Cross(y).Dot(z) < 0.5) { + LogDebug("magneticfield::volumeHandle") << "*** WARNING: Rotation is not RH "; + } + } + + // The global rotation + Surface::RotationType rotResult(float(x.X()), + float(x.Y()), + float(x.Z()), + float(y.X()), + float(y.Y()), + float(y.Z()), + float(z.X()), + float(z.Y()), + float(z.Z())); + + refPlane = new GloballyPositioned(posResult, rotResult); + + // Check correct orientation + if (debug) { + LogDebug("magneticfield::volumeHandle") << "Refplane pos " << refPlane->position(); + + // See comments above for the conventions for orientation. + LocalVector globalZdir(0., 0., 1.); // Local direction of the axis along global Z + + /* Preserve in case pseudotrap is needed again + if (shape_ == DDSolidShape::ddpseudotrap) { + globalZdir = LocalVector(0., 1., 0.); + } + */ + if (refPlane->toGlobal(globalZdir).z() < 0.) { + globalZdir = -globalZdir; + } + float chk = refPlane->toGlobal(globalZdir).dot(GlobalVector(0, 0, 1)); + if (chk < .999) + LogDebug("magneticfield::volumeHandle") << "*** WARNING RefPlane check failed!***" << chk; + } +} + +std::vector volumeHandle::sides() const { + std::vector result; + for (int i = 0; i < 6; ++i) { + // If this is just a master volume out of wich a 2pi volume + // should be built (e.g. central cylinder), skip the phi boundaries. + if (expand && (i == phiplus || i == phiminus)) + continue; + + // FIXME: Skip null inner degenerate cylindrical surface + if (shape_ == DDSolidShape::ddtubs && i == SurfaceOrientation::inner && theRMin < 0.001) + continue; + + ReferenceCountingPointer s = const_cast(surfaces[i].get()); + result.push_back(VolumeSide(s, GlobalFace(i), surfaces[i]->side(center_, 0.3))); + } + return result; +} + +// To allow the following code to work with both old and new DD, +// we need this function that merely returns its argument. +// For the old DD, another version of this function converts mm to cm. + +template +inline constexpr NumType convertUnits(NumType centimeters) { + return (centimeters); +} + +#include "MagneticField/GeomBuilder/src/buildBox.icc" +#include "MagneticField/GeomBuilder/src/buildTrap.icc" +#include "MagneticField/GeomBuilder/src/buildTubs.icc" +#include "MagneticField/GeomBuilder/src/buildCons.icc" +#include "MagneticField/GeomBuilder/src/buildTruncTubs.icc" diff --git a/MagneticField/GeomBuilder/plugins/dd4hep/volumeHandle.h b/MagneticField/GeomBuilder/plugins/dd4hep/volumeHandle.h new file mode 100644 index 0000000000000..4f90a74039f11 --- /dev/null +++ b/MagneticField/GeomBuilder/plugins/dd4hep/volumeHandle.h @@ -0,0 +1,60 @@ +#ifndef volumeHandle_H +#define volumeHandle_H + +/** \class volumeHandle + * A temporary container to cache info on a six-surface volume during + * the processing. Used to sort, organise, and build shared planes. + * One instance is created for each volume. The parameters of the + * boundary surfaces are calculated during construction. + * + * \author N. Amapane - INFN Torino + */ + +#include "DataFormats/GeometrySurface/interface/Surface.h" +#include "DetectorDescription/DDCMS/interface/DDFilteredView.h" +#include "DetectorDescription/DDCMS/interface/DDShapes.h" +#include "MagneticField/VolumeGeometry/interface/VolumeSide.h" +#include "MagneticField/GeomBuilder/src/BaseVolumeHandle.h" + +namespace magneticfield { + + typedef const char* ShapeType; + + class volumeHandle : public BaseVolumeHandle { + public: + volumeHandle(const cms::DDFilteredView& fv, bool expand2Pi = false, bool debugVal = false); + + // Disallow Default/copy ctor & assignment op. + // (we want to handle only pointers!!!) + volumeHandle(const volumeHandle& v) = delete; + volumeHandle operator=(const volumeHandle& v) = delete; + + // Shape at initialization + DDSolidShape shape() const override { return (shape_); } + + /// The surfaces and they orientation, as required to build a MagVolume. + std::vector sides() const override; + + private: + // initialise the refPlane + void referencePlane(const cms::DDFilteredView& fv); + + // Build the surfaces for a box + void buildBox(); + // Build the surfaces for a trapezoid + void buildTrap(); + // Build the surfaces for a ddtubs shape + void buildTubs(); + // Build the surfaces for a ddcons shape + void buildCons(); + // Build the surfaces for a ddtrunctubs shape + void buildTruncTubs(); + + // Shape at initialization + const DDSolidShape shape_; + const cms::DDFilteredView& solid; + // "solid" name is for backwards compatibility. Can be changed to "fview" after DD4hep migration. + }; +} // namespace magneticfield + +#endif diff --git a/MagneticField/GeomBuilder/src/BaseVolumeHandle.cc b/MagneticField/GeomBuilder/src/BaseVolumeHandle.cc new file mode 100644 index 0000000000000..aa02bc71e4807 --- /dev/null +++ b/MagneticField/GeomBuilder/src/BaseVolumeHandle.cc @@ -0,0 +1,221 @@ +/* + * See header file for a description of this class. + * + * \author N. Amapane - INFN Torino + */ + +#include "MagneticField/GeomBuilder/src/BaseVolumeHandle.h" + +#include "DataFormats/GeometrySurface/interface/Plane.h" +#include "DataFormats/GeometrySurface/interface/Cylinder.h" +#include "DataFormats/GeometrySurface/interface/Cone.h" +#include "DataFormats/GeometryVector/interface/CoordinateSets.h" + +#include "CLHEP/Units/GlobalSystemOfUnits.h" + +#include "MagneticField/Layers/interface/MagVerbosity.h" + +#include +#include +#include +#include +#include + +using namespace SurfaceOrientation; +using namespace std; +using namespace magneticfield; + +BaseVolumeHandle::BaseVolumeHandle(bool debugVal) + : magVolume(nullptr), + masterSector(1), + theRN(0.), + theRMin(0.), + theRMax(0.), + refPlane(nullptr), + isIronFlag(false), + debug(debugVal) {} + +BaseVolumeHandle::~BaseVolumeHandle() { + if (refPlane != nullptr) { + delete refPlane; + refPlane = nullptr; + } +} + +const Surface::GlobalPoint &BaseVolumeHandle::center() const { return center_; } + +void BaseVolumeHandle::buildPhiZSurf(double startPhi, double deltaPhi, double zhalf, double rCentr) { + // This is 100% equal for cons and tubs!!! + + GlobalVector planeXAxis = refPlane->toGlobal(LocalVector(1, 0, 0)); + GlobalVector planeYAxis = refPlane->toGlobal(LocalVector(0, 1, 0)); + GlobalVector planeZAxis = refPlane->toGlobal(LocalVector(0, 0, 1)); + + // Local Y axis of the faces at +-phi. + GlobalVector y_phiplus = refPlane->toGlobal(LocalVector(cos(startPhi + deltaPhi), sin(startPhi + deltaPhi), 0.)); + GlobalVector y_phiminus = refPlane->toGlobal(LocalVector(cos(startPhi), sin(startPhi), 0.)); + + Surface::RotationType rot_Z(planeXAxis, planeYAxis); + Surface::RotationType rot_phiplus(planeZAxis, y_phiplus); + Surface::RotationType rot_phiminus(planeZAxis, y_phiminus); + + GlobalPoint pos_zplus(center_.x(), center_.y(), center_.z() + zhalf); + GlobalPoint pos_zminus(center_.x(), center_.y(), center_.z() - zhalf); + // BEWARE: in this case, the origin for phiplus,phiminus surfaces is + // at radius R and not on a plane passing by center_ orthogonal to the radius. + GlobalPoint pos_phiplus( + refPlane->toGlobal(LocalPoint(rCentr * cos(startPhi + deltaPhi), rCentr * sin(startPhi + deltaPhi), 0.))); + GlobalPoint pos_phiminus(refPlane->toGlobal(LocalPoint(rCentr * cos(startPhi), rCentr * sin(startPhi), 0.))); + surfaces[zplus] = new Plane(pos_zplus, rot_Z); + surfaces[zminus] = new Plane(pos_zminus, rot_Z); + surfaces[phiplus] = new Plane(pos_phiplus, rot_phiplus); + surfaces[phiminus] = new Plane(pos_phiminus, rot_phiminus); + + if (debug) { + cout << "Actual Center at: " << center_ << " R " << center_.perp() << " phi " << center_.phi() << endl; + cout << "RN " << theRN << endl; + + cout << "pos_zplus " << pos_zplus << " " << pos_zplus.perp() << " " << pos_zplus.phi() << endl + << "pos_zminus " << pos_zminus << " " << pos_zminus.perp() << " " << pos_zminus.phi() << endl + << "pos_phiplus " << pos_phiplus << " " << pos_phiplus.perp() << " " << pos_phiplus.phi() << endl + << "pos_phiminus " << pos_phiminus << " " << pos_phiminus.perp() << " " << pos_phiminus.phi() << endl; + + cout << "y_phiplus " << y_phiplus << endl; + cout << "y_phiminus " << y_phiminus << endl; + + cout << "rot_Z " << surfaces[zplus]->toGlobal(LocalVector(0., 0., 1.)) << endl + << "rot_phi+ " << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " + << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl + << "rot_phi- " << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " + << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl; + } + + // // Check ordering. + if (debug) { + if (pos_zplus.z() < pos_zminus.z()) { + cout << "*** WARNING: pos_zplus < pos_zminus " << endl; + } + if (Geom::Phi(pos_phiplus.phi() - pos_phiminus.phi()) < 0.) { + cout << "*** WARNING: pos_phiplus < pos_phiminus " << endl; + } + } +} + +bool BaseVolumeHandle::sameSurface(const Surface &s1, Sides which_side, float tolerance) { + //Check for null comparison + if (&s1 == (surfaces[which_side]).get()) { + if (debug) + cout << " sameSurface: OK (same ptr)" << endl; + return true; + } + + const float maxtilt = 0.999; + + const Surface &s2 = *(surfaces[which_side]); + // Try with a plane. + const Plane *p1 = dynamic_cast(&s1); + if (p1 != nullptr) { + const Plane *p2 = dynamic_cast(&s2); + if (p2 == nullptr) { + if (debug) + cout << " sameSurface: different types" << endl; + return false; + } + + if ((fabs(p1->normalVector().dot(p2->normalVector())) > maxtilt) && + (fabs((p1->toLocal(p2->position())).z()) < tolerance)) { + if (debug) + cout << " sameSurface: OK " << fabs(p1->normalVector().dot(p2->normalVector())) << " " + << fabs((p1->toLocal(p2->position())).z()) << endl; + return true; + } else { + if (debug) + cout << " sameSurface: not the same: " << p1->normalVector() << p1->position() << endl + << " " << p2->normalVector() << p2->position() << endl + << fabs(p1->normalVector().dot(p2->normalVector())) << " " << (p1->toLocal(p2->position())).z() << endl; + return false; + } + } + + // Try with a cylinder. + const Cylinder *cy1 = dynamic_cast(&s1); + if (cy1 != nullptr) { + const Cylinder *cy2 = dynamic_cast(&s2); + if (cy2 == nullptr) { + if (debug) + cout << " sameSurface: different types" << endl; + return false; + } + // Assume axis is the same! + if (fabs(cy1->radius() - cy2->radius()) < tolerance) { + return true; + } else { + return false; + } + } + + // Try with a cone. + const Cone *co1 = dynamic_cast(&s1); + if (co1 != nullptr) { + const Cone *co2 = dynamic_cast(&s2); + if (co2 == nullptr) { + if (debug) + cout << " sameSurface: different types" << endl; + return false; + } + // FIXME + if (fabs(co1->openingAngle() - co2->openingAngle()) < maxtilt && + (co1->vertex() - co2->vertex()).mag() < tolerance) { + return true; + } else { + return false; + } + } + + if (debug) + cout << " sameSurface: unknown surfaces..." << endl; + return false; +} + +bool BaseVolumeHandle::setSurface(const Surface &s1, Sides which_side) { + //Check for null assignment + if (&s1 == (surfaces[which_side]).get()) { + isAssigned[which_side] = true; + return true; + } + + if (!sameSurface(s1, which_side)) { + cout << "***ERROR: setSurface: trying to assign a surface that does not match destination surface. Skipping." + << endl; + const Surface &s2 = *(surfaces[which_side]); + //FIXME: Just planes for the time being!!! + const Plane *p1 = dynamic_cast(&s1); + const Plane *p2 = dynamic_cast(&s2); + if (p1 != nullptr && p2 != nullptr) + cout << p1->normalVector() << p1->position() << endl << p2->normalVector() << p2->position() << endl; + return false; + } + + if (isAssigned[which_side]) { + if (&s1 != (surfaces[which_side]).get()) { + cout << "*** WARNING BaseVolumeHandle::setSurface: trying to reassign a surface to a different surface instance" + << endl; + return false; + } + } else { + surfaces[which_side] = &s1; + isAssigned[which_side] = true; + if (debug) + cout << " Volume " << name << " # " << copyno << " Assigned: " << (int)which_side << endl; + return true; + } + + return false; // let the compiler be happy +} + +const Surface &BaseVolumeHandle::surface(Sides which_side) const { return *(surfaces[which_side]); } + +const Surface &BaseVolumeHandle::surface(int which_side) const { + assert(which_side >= 0 && which_side < 6); + return *(surfaces[which_side]); +} diff --git a/MagneticField/GeomBuilder/src/BaseVolumeHandle.h b/MagneticField/GeomBuilder/src/BaseVolumeHandle.h new file mode 100644 index 0000000000000..65ea9693f5dcd --- /dev/null +++ b/MagneticField/GeomBuilder/src/BaseVolumeHandle.h @@ -0,0 +1,219 @@ +#ifndef MagneticField_GeomBuilder_BaseVolumeHandle_H +#define MagneticField_GeomBuilder_BaseVolumeHandle_H + +/** \class BaseVolumeHandle + * A temporary container to cache info on a six-surface volume during + * the processing. Used to sort, organise, and build shared planes. + * One instance is created for each DDVolume. The parameters of the + * boundary surfaces are calculated during construction. + * + * \author N. Amapane - INFN Torino + */ + +#include "DataFormats/GeometrySurface/interface/Surface.h" +#include "DataFormats/GeometrySurface/interface/ReferenceCounted.h" +#include "DetectorDescription/Core/interface/DDSolidShapes.h" +#include "MagneticField/VolumeGeometry/interface/VolumeSide.h" + +class MagVolume6Faces; + +namespace magneticfield { + + class BaseVolumeHandle { + public: + typedef Surface::GlobalPoint GlobalPoint; + typedef Surface::LocalPoint LocalPoint; + typedef Surface::LocalVector LocalVector; + typedef SurfaceOrientation::GlobalFace Sides; + + BaseVolumeHandle(bool debugVal = false); + + // Disallow Default/copy ctor + // (we want to handle only pointers!!!) + BaseVolumeHandle(const BaseVolumeHandle& v) = delete; + + virtual ~BaseVolumeHandle(); + + /// Return the center of the volume + const GlobalPoint& center() const; + + /// Distance of (x,y) plane from origin + const double RN() const { return theRN; } + + /// Get the current surface on specified side. + const Surface& surface(int which_side) const; + + const Surface& surface(Sides which_side) const; + + /// Find out if two surfaces are the same physical surface + bool sameSurface(const Surface& s1, Sides which_side, float tolerance = 0.01); + + /// Assign a shared surface perorming sanity checks. + bool setSurface(const Surface& s1, Sides which_side); + + /// if the specified surface has been matched. + bool isPlaneMatched(int which_side) const { return isAssigned[which_side]; } + + int references(int which_side) const { // FIXME! + /* return surfaces[which_side]->references(); */ + return 0; + } + + /// Name of the volume + std::string name; + + /// Name of magnetic field table file + std::string magFile; + + /// volume number + unsigned short volumeno; + + /// copy number + unsigned short copyno; + + // Phi ranges: Used by: LessDPhiMax; bSector; bSlab::[min|max]Phi(); + // MagBSector, MagBRod + + /// Minimum value of phi covered by the volume + // FIXME: actually returns phi of the point on median plane of the -phi + // surface, except for trapezoids where the absoulte min has been implemented + Geom::Phi minPhi() const { return thePhiMin; } + + /// Maximum value of phi covered by the volume + // FIXME: actually returns phi of the point on median plane of the +phi + // surface + Geom::Phi maxPhi() const { return surface(SurfaceOrientation::phiplus).position().phi(); } + + /// Z limits. + // ASSUMPTION: Computed on median Z plane, but this is not a problem since + // all Z planes are orthogonal to the beam line in the current geometry. + double minZ() const { return surface(SurfaceOrientation::zminus).position().z(); } + double maxZ() const { return surface(SurfaceOrientation::zplus).position().z(); } + + /// Minimum R for any point within the volume + double minR() const { return theRMin; } + + /// Position and rotation + const GloballyPositioned* placement() const { return refPlane; } + + /// The surfaces and they orientation, as required to build a MagVolume. + virtual std::vector sides() const = 0; + + /// Pointer to the final MagVolume (must be set from outside) + MagVolume6Faces* magVolume; + + bool toExpand() const { return expand; } + + /// Temporary hack to pass information on material. Will eventually be replaced! + bool isIron() const { return isIronFlag; } + + /// The sector for which an interpolator for this class of volumes should be built + int masterSector; + + /// Shape of the solid + virtual DDSolidShape shape() const = 0; + + protected: + typedef ConstReferenceCountingPointer RCPS; + + // The volume's six surfaces. + RCPS surfaces[6]; + // If the corresponding surface has been assigned to a shared surface. + bool isAssigned[6]; + + // Build phi, z surfaces (common for ddtubs and ddcons) + void buildPhiZSurf(double startPhi, double deltaPhi, double zhalf, double rCentr); + + // Distance from the origin along the normal to the volume's zphi plane. + double theRN; + + // Max and min radius for _any_ point within the volume + // FIXME! + double theRMin; + double theRMax; + Geom::Phi thePhiMin; + + // The refPlane is the "main plane" for the solid. It corresponds to the + // x,y plane in the DDD local frame, and defines a frame where the local + // coordinates are the same as in DDD. + GloballyPositioned* refPlane; + + // the center of the volume + GlobalPoint center_; + + // Flag this as a master volume out of wich a 2pi volume should be built + // (e.g. central cylinder); this is taken into account by sides(). + bool expand; + + // Temporary hack to keep information on material. Will eventually be replaced! + bool isIronFlag; + + const bool debug; + }; + + typedef std::vector handles; + + // Extractors for precomputed_value_sort() (safe sorting) + + // To sort volumes in Z + struct ExtractZ { + double operator()(const BaseVolumeHandle* v) const { return v->center().z(); } + }; + + // To sort volumes in abs(Z) + struct ExtractAbsZ { + double operator()(const BaseVolumeHandle* v) const { return fabs(v->center().z()); } + }; + + // To sort volumes in phi (from -pi to pi). + struct ExtractPhi { + double operator()(const BaseVolumeHandle* v) const { + // note that Geom::Phi is implicitly converted to double. + // Periodicity is guaranteed. + return v->center().phi(); + } + }; + + // To sort volumes based on max phi(from -pi to pi). + struct ExtractPhiMax { + double operator()(const BaseVolumeHandle* v) const { + // note that Geom::Phi is implicitly converted to double. + // Periodicity is guaranteed. + return v->maxPhi(); + } + }; + + // To sort volumes in R + struct ExtractR { + double operator()(const BaseVolumeHandle* v) const { return v->center().perp(); } + }; + + // To sort volumes in RN (distance of (x,y) plane from origin) + struct ExtractRN { + double operator()(const BaseVolumeHandle* v) const { return v->RN(); } + }; + + // To sort angles within any range SMALLER THAN PI "counter-clockwise", + // even if the angles cross the pi boundary. + // CAVEAT: // The result is undefined if the input values cover a + // range larger than pi!!! + struct LessDPhi { + bool operator()(double phi1, double phi2) const { + // handle periodicity + return ((Geom::Phi(phi2) - Geom::Phi(phi1)) > 0.); + } + }; + + // Compare the Z of volumes. + // Should be used ONLY for std::max_element and std::min_element + // and NEVER for sorting (use precomputed_value_sort with ExtractZ instead) + struct LessZ { + bool operator()(const BaseVolumeHandle* v1, const BaseVolumeHandle* v2) const { + if (v1->center().z() < v2->center().z()) + return true; + return false; + } + }; +} // namespace magneticfield + +#endif diff --git a/MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.cc b/MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.cc index fd80822e0e104..c49d143de0f4d 100644 --- a/MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.cc +++ b/MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.cc @@ -24,7 +24,6 @@ #include "DetectorDescription/Core/interface/DDFilter.h" #include "Utilities/BinningTools/interface/ClusterizingHistogram.h" -#include "CLHEP/Units/GlobalSystemOfUnits.h" #include "MagneticField/Interpolation/interface/MagProviderInterpol.h" #include "MagneticField/Interpolation/interface/MFGridFactory.h" @@ -35,6 +34,7 @@ #include "MagneticField/Layers/interface/MagVerbosity.h" #include "DataFormats/GeometryVector/interface/Pi.h" +#include "DataFormats/Math/interface/GeantUnits.h" #include "FWCore/Utilities/interface/Exception.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" @@ -51,14 +51,11 @@ #include #include "Utilities/General/interface/precomputed_value_sort.h" -bool MagGeoBuilderFromDDD::debug; - using namespace std; using namespace magneticfield; MagGeoBuilderFromDDD::MagGeoBuilderFromDDD(string tableSet_, int geometryVersion_, bool debug_) - : tableSet(tableSet_), geometryVersion(geometryVersion_), theGridFiles(nullptr) { - debug = debug_; + : tableSet(tableSet_), geometryVersion(geometryVersion_), theGridFiles(nullptr), debug(debug_) { if (debug) cout << "Constructing a MagGeoBuilderFromDDD" << endl; } @@ -315,14 +312,14 @@ void MagGeoBuilderFromDDD::build(const DDCompactView& cpva) { while ((*separ)->RN() < rSepar) ++separ; - bLayer thislayer(ringStart, separ); + bLayer thislayer(ringStart, separ, debug); layers.push_back(thislayer); ringStart = separ; } { if (debug) cout << " Layer at RN = " << rClust.back(); - bLayer thislayer(separ, last); + bLayer thislayer(separ, last, debug); layers.push_back(thislayer); } @@ -378,7 +375,7 @@ void MagGeoBuilderFromDDD::build(const DDCompactView& cpva) { } } - sectors.push_back(eSector(eVolumes.begin() + ((i)*offset), eVolumes.begin() + ((i + 1) * offset))); + sectors.push_back(eSector(eVolumes.begin() + ((i)*offset), eVolumes.begin() + ((i + 1) * offset), debug)); } if (debug) diff --git a/MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.h b/MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.h index f665dec6dcdc4..e1564b008e131 100644 --- a/MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.h +++ b/MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.h @@ -8,7 +8,6 @@ * * \author N. Amapane - INFN Torino */ -#include "DataFormats/GeometrySurface/interface/ReferenceCounted.h" #include "MagneticField/Interpolation/interface/MagProviderInterpol.h" #include "DetectorDescription/Core/interface/DDCompactView.h" @@ -27,6 +26,8 @@ namespace magneticfield { class VolumeBasedMagneticFieldESProducer; class VolumeBasedMagneticFieldESProducerFromDB; class AutoMagneticFieldESProducer; + class BaseVolumeHandle; // Needs to be public to share code with DD4hep + typedef std::vector handles; } // namespace magneticfield class MagGeoBuilderFromDDD { @@ -54,9 +55,10 @@ class MagGeoBuilderFromDDD { float maxZ() const; -private: - typedef ConstReferenceCountingPointer RCPS; + // Temporary container to manipulate volumes and their surfaces. + class volumeHandle; // Needs to be public to share code with DD4hep +private: // Build the geometry. //virtual void build(); virtual void build(const DDCompactView& cpv); @@ -71,48 +73,21 @@ class MagGeoBuilderFromDDD { std::vector barrelVolumes() const; std::vector endcapVolumes() const; - // Temporary container to manipulate volumes and their surfaces. - class volumeHandle; - typedef std::vector handles; - // Build interpolator for the volume with "correct" rotation void buildInterpolator(const volumeHandle* vol, std::map& interpolators); // Build all MagVolumes setting the MagProviderInterpol - void buildMagVolumes(const handles& volumes, std::map& interpolators); + void buildMagVolumes(const magneticfield::handles& volumes, + std::map& interpolators); // Print checksums for surfaces. - void summary(handles& volumes); + void summary(magneticfield::handles& volumes); // Perform simple sanity checks - void testInside(handles& volumes); - - // A layer of barrel volumes. - class bLayer; - // A sector of volumes in a layer. - class bSector; - // A rod of volumes in a sector. - class bRod; - // A slab of volumes in a rod. - class bSlab; - // A sector of endcap volumes. - class eSector; - // A layer of volumes in an endcap sector. - class eLayer; - - // Extractors for precomputed_value_sort (to sort containers of volumeHandles). - struct ExtractZ; - struct ExtractAbsZ; - struct ExtractPhi; - struct ExtractPhiMax; - struct ExtractR; - struct ExtractRN; - struct LessDPhi; - // This one to be used only for max_element and min_element - struct LessZ; - - handles bVolumes; // the barrel volumes. - handles eVolumes; // the endcap volumes. + void testInside(magneticfield::handles& volumes); + + magneticfield::handles bVolumes; // the barrel volumes. + magneticfield::handles eVolumes; // the endcap volumes. std::vector mBLayers; // Finally built barrel geometry std::vector mESectors; // Finally built barrel geometry @@ -123,6 +98,6 @@ class MagGeoBuilderFromDDD { std::map theScalingFactors; const magneticfield::TableFileMap* theGridFiles; // Non-owned pointer assumed to be valid until build() is called - static bool debug; + const bool debug; }; #endif diff --git a/MagneticField/GeomBuilder/src/bLayer.cc b/MagneticField/GeomBuilder/src/bLayer.cc index 23ea5b95e9fb9..fadaf39d47427 100644 --- a/MagneticField/GeomBuilder/src/bLayer.cc +++ b/MagneticField/GeomBuilder/src/bLayer.cc @@ -7,6 +7,7 @@ */ #include "MagneticField/GeomBuilder/src/bLayer.h" +#include "MagneticField/GeomBuilder/src/printUniqueNames.h" #include "MagneticField/VolumeGeometry/interface/MagVolume6Faces.h" #include "MagneticField/Layers/interface/MagBLayer.h" #include "MagneticField/Layers/interface/MagVerbosity.h" @@ -14,21 +15,22 @@ #include "Utilities/General/interface/precomputed_value_sort.h" using namespace SurfaceOrientation; +using namespace magneticfield; //The ctor is in charge of finding sectors inside the layer. -MagGeoBuilderFromDDD::bLayer::bLayer(handles::const_iterator begin, handles::const_iterator end) - : size(end - begin), theVolumes(begin, end), mlayer(nullptr) { +bLayer::bLayer(handles::const_iterator begin, handles::const_iterator end, bool debugFlag) + : size(end - begin), theVolumes(begin, end), mlayer(nullptr), debug(debugFlag) { // Sort in phi precomputed_value_sort(theVolumes.begin(), theVolumes.end(), ExtractPhi()); - if (MagGeoBuilderFromDDD::debug) { + if (debug) { std::cout << " elements: " << theVolumes.size() << " unique volumes: "; - volumeHandle::printUniqueNames(theVolumes.begin(), theVolumes.end()); + printUniqueNames(theVolumes.begin(), theVolumes.end()); } // Find sectors in phi handles::iterator secBegin = theVolumes.begin(); - handles::iterator secEnd; + handles::iterator secEnd = secBegin; int binOffset = 0; const Surface& refSurf = (*secBegin)->surface(outer); @@ -45,7 +47,7 @@ MagGeoBuilderFromDDD::bLayer::bLayer(handles::const_iterator begin, handles::con if (size == 1) { // Only one volume; this is the case for barrel // cylinders. // FIXME sectors.push_back(bSector(theVolumes.begin(),theVolumes.end()); - if (MagGeoBuilderFromDDD::debug) + if (debug) std::cout << " Sector is just one volume." << std::endl; } else if (size == 12 || // In this case, each volume is a sector. @@ -55,7 +57,7 @@ MagGeoBuilderFromDDD::bLayer::bLayer(handles::const_iterator begin, handles::con } else { // there are more than one volume per sector. float tolerance = 0.025; // 250 micron do { - if (MagGeoBuilderFromDDD::debug) + if (debug) std::cout << (*secBegin)->name << " " << (*secBegin)->copyno << std::endl; ++secBegin; } while ( @@ -83,7 +85,7 @@ MagGeoBuilderFromDDD::bLayer::bLayer(handles::const_iterator begin, handles::con } } - if (MagGeoBuilderFromDDD::debug) { + if (debug) { std::cout << " First sector: volumes " << secEnd - theVolumes.begin() << " from " << newbegin << " (phi = " << (*secBegin)->center().phi() << ") " << " to " << newend << " (phi = " << (*secEnd)->center().phi() << ") " @@ -97,42 +99,40 @@ MagGeoBuilderFromDDD::bLayer::bLayer(handles::const_iterator begin, handles::con for (int i = 0; i < 12; ++i) { int isec = (i + binOffset) % 12; sectors[isec >= 0 ? isec : isec + 12] = - bSector(theVolumes.begin() + ((i)*offset), theVolumes.begin() + ((i + 1) * offset)); + bSector(theVolumes.begin() + ((i)*offset), theVolumes.begin() + ((i + 1) * offset), debug); } } - if (MagGeoBuilderFromDDD::debug) + if (debug) std::cout << "-----------------------" << std::endl; } -MagGeoBuilderFromDDD::bLayer::~bLayer() {} - -int MagGeoBuilderFromDDD::bLayer::bin(int i) const { +int bLayer::bin(int i) const { i = i % size; return (i >= 0 ? i : i + size); } -// const MagGeoBuilderFromDDD::bSector & -// MagGeoBuilderFromDDD::bLayer::sector(int i) const { +// const bSector & +// bLayer::sector(int i) const { // i = i%12; // return sectors[i>=0?i:i+12]; // } -double MagGeoBuilderFromDDD::bLayer::minR() const { +double bLayer::minR() const { // ASSUMPTION: a layer is only 1 volume thick (by construction). return theVolumes.front()->minR(); } -// double MagGeoBuilderFromDDD::bLayer::maxR() const { +// double bLayer::maxR() const { // // ASSUMPTION: a layer is only 1 volume thick (by construction). // return theVolumes.front()->maxR(); // } -MagBLayer* MagGeoBuilderFromDDD::bLayer::buildMagBLayer() const { +MagBLayer* bLayer::buildMagBLayer() const { if (mlayer == nullptr) { // If we have only one volume, do not build any MagBSector. if (sectors.empty()) { - if (MagGeoBuilderFromDDD::debug && size != 0) { + if (debug && size != 0) { std::cout << "ERROR: bLayer::buildMagBLayer, 0 sectors but " << size << " volumes" << std::endl; } // Technically we might have only one bSector built and we would diff --git a/MagneticField/GeomBuilder/src/bLayer.h b/MagneticField/GeomBuilder/src/bLayer.h index 0eb19446b9c5a..615910786abef 100644 --- a/MagneticField/GeomBuilder/src/bLayer.h +++ b/MagneticField/GeomBuilder/src/bLayer.h @@ -1,55 +1,58 @@ #ifndef bLayer_H #define bLayer_H -/** \class MagGeoBuilderFromDDD::bLayer +/** \class bLayer * A layer of barrel volumes. Holds a list of volumes and 12 sectors. * It is assumed that the geometry is 12-fold periodic in phi! * * \author N. Amapane - INFN Torino */ -#include "MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.h" -#include "MagneticField/GeomBuilder/src/volumeHandle.h" #include "MagneticField/GeomBuilder/src/bSector.h" class MagBLayer; -class MagGeoBuilderFromDDD::bLayer { -public: - /// Constructor from list of volumes - bLayer(handles::const_iterator begin, handles::const_iterator end); +namespace magneticfield { - /// Destructor - ~bLayer(); + class bLayer { + public: + /// Constructor from list of volumes + bLayer(handles::const_iterator begin, handles::const_iterator end, bool debugFlag = false); - /// Distance from center along normal of sectors. - const float RN() const { return theVolumes.front()->RN(); } + /// Destructor + ~bLayer() = default; - /// Return the list of all volumes. - const handles& volumes() const { return theVolumes; } + /// Distance from center along normal of sectors. + const float RN() const { return theVolumes.front()->RN(); } - /// Return sector at i (handling periodicity) - // const bSector & sector(int i) const; + /// Return the list of all volumes. + const handles& volumes() const { return theVolumes; } - /// Min R (conservative guess). - double minR() const; + /// Return sector at i (handling periodicity) + // const bSector & sector(int i) const; - // Depends on volumeHandle::maxR(), which actually returns max RN. - // (should be changed?) - // double maxR() const; + /// Min R (conservative guess). + double minR() const; - /// Construct the MagBLayer upon request. - MagBLayer* buildMagBLayer() const; + // Depends on volumeHandle::maxR(), which actually returns max RN. + // (should be changed?) + // double maxR() const; -private: - int size; //< the number of volumes + /// Construct the MagBLayer upon request. + MagBLayer* buildMagBLayer() const; - // Check periodicity; - int bin(int i) const; + private: + int size; //< the number of volumes - std::vector sectors; // the sectors in this layer - handles theVolumes; // pointer to all volumes in this layer + // Check periodicity; + int bin(int i) const; + + std::vector sectors; // the sectors in this layer + handles theVolumes; // pointer to all volumes in this layer + + mutable MagBLayer* mlayer; + const bool debug; + }; +} // namespace magneticfield - mutable MagBLayer* mlayer; -}; #endif diff --git a/MagneticField/GeomBuilder/src/bRod.cc b/MagneticField/GeomBuilder/src/bRod.cc index eb3d2431652bb..cd8615180d39a 100644 --- a/MagneticField/GeomBuilder/src/bRod.cc +++ b/MagneticField/GeomBuilder/src/bRod.cc @@ -7,18 +7,18 @@ */ #include "MagneticField/GeomBuilder/src/bRod.h" +#include "MagneticField/GeomBuilder/src/printUniqueNames.h" #include "Utilities/BinningTools/interface/ClusterizingHistogram.h" #include "MagneticField/Layers/interface/MagBRod.h" #include "MagneticField/Layers/interface/MagVerbosity.h" #include "Utilities/General/interface/precomputed_value_sort.h" using namespace SurfaceOrientation; - -MagGeoBuilderFromDDD::bRod::~bRod() {} +using namespace magneticfield; //The ctor is in charge of finding slabs inside the rod. -MagGeoBuilderFromDDD::bRod::bRod(handles::const_iterator begin, handles::const_iterator end) - : volumes(begin, end), mrod(nullptr) { +bRod::bRod(handles::const_iterator begin, handles::const_iterator end, bool debugVal) + : volumes(begin, end), mrod(nullptr), debug(debugVal) { precomputed_value_sort(volumes.begin(), volumes.end(), ExtractZ()); // Clusterize in Z @@ -27,7 +27,7 @@ MagGeoBuilderFromDDD::bRod::bRod(handles::const_iterator begin, handles::const_i float zmax = volumes.back()->center().z() + resolution; ClusterizingHistogram hisZ(int((zmax - zmin) / resolution) + 1, zmin, zmax); - if (MagGeoBuilderFromDDD::debug) + if (debug) std::cout << " Z slabs: " << zmin << " " << zmax << std::endl; handles::const_iterator first = volumes.begin(); @@ -38,7 +38,7 @@ MagGeoBuilderFromDDD::bRod::bRod(handles::const_iterator begin, handles::const_i } std::vector zClust = hisZ.clusterize(resolution); - if (MagGeoBuilderFromDDD::debug) + if (debug) std::cout << " Found " << zClust.size() << " clusters in Z, " << " slabs: " << std::endl; @@ -49,20 +49,20 @@ MagGeoBuilderFromDDD::bRod::bRod(handles::const_iterator begin, handles::const_i float zSepar = (zClust[i] + zClust[i + 1]) / 2.f; while ((*separ)->center().z() < zSepar) ++separ; - if (MagGeoBuilderFromDDD::debug) { + if (debug) { std::cout << " Slab at: " << zClust[i] << " elements: " << separ - slabStart << " unique volumes: "; - volumeHandle::printUniqueNames(slabStart, separ); + printUniqueNames(slabStart, separ); } - slabs.push_back(bSlab(slabStart, separ)); + slabs.push_back(bSlab(slabStart, separ, debug)); slabStart = separ; } { - if (MagGeoBuilderFromDDD::debug) { + if (debug) { std::cout << " Slab at: " << zClust.back() << " elements: " << last - separ << " unique volumes: "; - volumeHandle::printUniqueNames(separ, last); + printUniqueNames(separ, last); } - slabs.push_back(bSlab(separ, last)); + slabs.push_back(bSlab(separ, last, debug)); } // Check that all slabs have the same dphi. @@ -71,13 +71,15 @@ MagGeoBuilderFromDDD::bRod::bRod(handles::const_iterator begin, handles::const_i Geom::Phi phimin = (*i).minPhi(); for (++i; i != slabs.end(); ++i) { if (fabs(phimax - (*i).maxPhi()) > 0.001 || fabs(phimin - (*i).minPhi()) > 0.001) { - if (MagGeoBuilderFromDDD::debug) - std::cout << "*** WARNING: slabs in this rod have different dphi!" << std::endl; + if (debug) { + std::cout << "*** WARNING: slabs in this rod have different dphi! minphi " << phimin; + std::cout << " != " << (*i).minPhi() << " or maxphi " << phimax << " != " << (*i).maxPhi() << std::endl; + } } } } -MagBRod* MagGeoBuilderFromDDD::bRod::buildMagBRod() const { +MagBRod* bRod::buildMagBRod() const { if (mrod == nullptr) { std::vector mSlabs; for (std::vector::const_iterator slab = slabs.begin(); slab != slabs.end(); ++slab) { diff --git a/MagneticField/GeomBuilder/src/bRod.h b/MagneticField/GeomBuilder/src/bRod.h index 93a1db362d42e..05740f6457b7d 100644 --- a/MagneticField/GeomBuilder/src/bRod.h +++ b/MagneticField/GeomBuilder/src/bRod.h @@ -8,30 +8,32 @@ * \author N. Amapane - INFN Torino */ -#include "MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.h" -#include "MagneticField/GeomBuilder/src/volumeHandle.h" #include "MagneticField/GeomBuilder/src/bSlab.h" class MagBRod; -class MagGeoBuilderFromDDD::bRod { -public: - /// Constructor from list of volumes - bRod(handles::const_iterator begin, handles::const_iterator end); +namespace magneticfield { - /// Destructor - ~bRod(); + class bRod { + public: + /// Constructor from list of volumes + bRod(handles::const_iterator begin, handles::const_iterator end, bool debugVal = false); - /// Distance from center along sector normal. - const float RN() const { return volumes.front()->RN(); } + /// Destructor + ~bRod() = default; - /// Construct the MagBRod upon request. - MagBRod* buildMagBRod() const; + /// Distance from center along sector normal. + const float RN() const { return volumes.front()->RN(); } -private: - std::vector slabs; - handles volumes; // pointers to all volumes in the rod - mutable MagBRod* mrod; -}; + /// Construct the MagBRod upon request. + MagBRod* buildMagBRod() const; + + private: + std::vector slabs; + handles volumes; // pointers to all volumes in the rod + mutable MagBRod* mrod; + bool debug; // Allow assignment from other bRod objects + }; +} // namespace magneticfield #endif diff --git a/MagneticField/GeomBuilder/src/bSector.cc b/MagneticField/GeomBuilder/src/bSector.cc index eaece6b43594d..deeddf3515986 100644 --- a/MagneticField/GeomBuilder/src/bSector.cc +++ b/MagneticField/GeomBuilder/src/bSector.cc @@ -6,35 +6,35 @@ * \author N. Amapane - INFN Torino */ +#include "DataFormats/GeometrySurface/interface/Surface.h" #include "MagneticField/GeomBuilder/src/bSector.h" #include "Utilities/BinningTools/interface/ClusterizingHistogram.h" #include "MagneticField/Layers/interface/MagBSector.h" #include "MagneticField/Layers/interface/MagVerbosity.h" - +#include "MagneticField/GeomBuilder/src/printUniqueNames.h" #include "Utilities/General/interface/precomputed_value_sort.h" #include using namespace SurfaceOrientation; using namespace std; +using namespace magneticfield; // Default ctor needed to have arrays. -MagGeoBuilderFromDDD::bSector::bSector() {} - -MagGeoBuilderFromDDD::bSector::~bSector() {} +bSector::bSector() : debug(false) {} // The ctor is in charge of finding rods inside the sector. -MagGeoBuilderFromDDD::bSector::bSector(handles::const_iterator begin, handles::const_iterator end) - : volumes(begin, end), msector(nullptr) { - if (MagGeoBuilderFromDDD::debug) +bSector::bSector(handles::const_iterator begin, handles::const_iterator end, bool debugVal) + : volumes(begin, end), msector(nullptr), debug(debugVal) { + if (debug) cout << " Sector at Phi " << volumes.front()->center().phi() << " " << volumes.back()->center().phi() << endl; if (volumes.size() == 1) { - if (MagGeoBuilderFromDDD::debug) { + if (debug) { cout << " Rod at: 0 elements: " << end - begin << " unique volumes: "; - volumeHandle::printUniqueNames(begin, end); + printUniqueNames(begin, end); } - rods.push_back(bRod(begin, end)); + rods.push_back(bRod(begin, end, debug)); } else { // Clusterize in phi. Use bin edge so that complete clusters can be // easily found (not trivial using bin centers!) @@ -46,11 +46,21 @@ MagGeoBuilderFromDDD::bSector::bSector(handles::const_iterator begin, handles::c // Sort volumes in DELTA phi - i.e. phi(j)-phi(i) > 0 if j>1. precomputed_value_sort(volumes.begin(), volumes.end(), ExtractPhiMax(), LessDPhi()); - const Geom::Phi resolution(0.01); // rad + const float resolution(0.01); // rad Geom::Phi phi0 = volumes.front()->maxPhi(); - float phiMin = -(float)resolution; - float phiMax = volumes.back()->maxPhi() - phi0 + resolution; ///FIXME: (float) resolution; ?? - + float phiMin = -resolution; ///FIXME: (float) resolution; ?? + + // Careful -- Phi overloads arithmetic operators and will wrap around each step of a calculation, + // so use casts to prevent intermediate value wrap-arounds. + float phiTmp = static_cast(volumes.back()->maxPhi()) - static_cast(phi0) + resolution; + const float phiMax = angle0to2pi::make0To2pi(phiTmp); // Ensure 0-2pi + + if (debug) { + cout << "volumes size = " << volumes.size(); + cout << ", phi0 = " << phi0 << ", volumes.back()->maxPhi() = " << volumes.back()->maxPhi(); + cout << ", phiMin = " << phiMin << ", phiMax = " << phiMax; + cout << ", int((phiMax - phiMin) / resolution) + 1 = " << int((phiMax - phiMin) / resolution) + 1 << endl; + } ClusterizingHistogram hisPhi(int((phiMax - phiMin) / resolution) + 1, phiMin, phiMax); handles::const_iterator first = volumes.begin(); @@ -61,7 +71,7 @@ MagGeoBuilderFromDDD::bSector::bSector(handles::const_iterator begin, handles::c } vector phiClust = hisPhi.clusterize(resolution); - if (MagGeoBuilderFromDDD::debug) + if (debug) cout << " Found " << phiClust.size() << " clusters in Phi, " << " rods: " << endl; @@ -78,11 +88,11 @@ MagGeoBuilderFromDDD::bSector::bSector(handles::const_iterator begin, handles::c } else { phiSepar = phiMax; } - if (MagGeoBuilderFromDDD::debug) + if (debug) cout << " cluster " << i << " phisepar " << phiSepar << endl; while (separ < last && (*separ)->maxPhi() - phi0 < phiSepar) { DZ1 += ((*separ)->maxZ() - (*separ)->minZ()); - if (MagGeoBuilderFromDDD::debug) + if (debug) cout << " " << (*separ)->name << " " << (*separ)->maxPhi() - phi0 << " " << (*separ)->maxZ() << " " << (*separ)->minZ() << " " << DZ1 << endl; ++separ; @@ -91,25 +101,25 @@ MagGeoBuilderFromDDD::bSector::bSector(handles::const_iterator begin, handles::c // FIXME: print warning for small discrepancies. Tolerance (below) // had to be increased since discrepancies sum to up to ~ 2 mm. if (fabs(DZ - DZ1) > 0.001 && fabs(DZ - DZ1) < 0.5) { - if (MagGeoBuilderFromDDD::debug) + if (debug) cout << "*** WARNING: Z lenght mismatch by " << DZ - DZ1 << " " << DZ << " " << DZ1 << endl; } if (fabs(DZ - DZ1) > 0.25) { // FIXME hardcoded tolerance - if (MagGeoBuilderFromDDD::debug) + if (debug) cout << " Incomplete, use also next cluster: " << DZ << " " << DZ1 << " " << DZ - DZ1 << endl; DZ1 = 0.; continue; } else if (DZ1 > DZ + 0.05) { // Wrong: went past max lenght // FIXME hardcoded tolerance cout << " *** ERROR: bSector finding messed up." << endl; - volumeHandle::printUniqueNames(rodStart, separ); + printUniqueNames(rodStart, separ); DZ1 = 0.; } else { - if (MagGeoBuilderFromDDD::debug) { + if (debug) { cout << " Rod at: " << phiClust[i] << " elements: " << separ - rodStart << " unique volumes: "; - volumeHandle::printUniqueNames(rodStart, separ); + printUniqueNames(rodStart, separ); } - rods.push_back(bRod(rodStart, separ)); + rods.push_back(bRod(rodStart, separ, debug)); rodStart = separ; DZ1 = 0.; } @@ -117,18 +127,19 @@ MagGeoBuilderFromDDD::bSector::bSector(handles::const_iterator begin, handles::c if (rods.empty()) cout << " *** ERROR: bSector has no rods " << DZ << " " << DZ1 << endl; - if (MagGeoBuilderFromDDD::debug) + if (debug) cout << "-----------------------" << endl; } } -MagBSector* MagGeoBuilderFromDDD::bSector::buildMagBSector() const { +MagBSector* bSector::buildMagBSector() const { if (msector == nullptr) { vector mRods; for (vector::const_iterator rod = rods.begin(); rod != rods.end(); ++rod) { mRods.push_back((*rod).buildMagBRod()); } msector = new MagBSector(mRods, volumes.front()->minPhi()); //FIXME + // Never deleted. When is it safe to delete it? } return msector; } diff --git a/MagneticField/GeomBuilder/src/bSector.h b/MagneticField/GeomBuilder/src/bSector.h index 1724bc1d96941..ff690ec5cf0b3 100644 --- a/MagneticField/GeomBuilder/src/bSector.h +++ b/MagneticField/GeomBuilder/src/bSector.h @@ -1,42 +1,45 @@ #ifndef bSector_H #define bSector_H -/** \class MagGeoBuilderFromDDD::bSector +/** \class bSector * A sector of volumes in a barrel layer (i.e. only 1 element in R) * One sector is composed of 1 or more rods. * * \author N. Amapane - INFN Torino */ -#include "MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.h" -#include "MagneticField/GeomBuilder/src/volumeHandle.h" #include "MagneticField/GeomBuilder/src/bRod.h" class MagBSector; -class MagGeoBuilderFromDDD::bSector { -public: - /// Default ctor is needed to have arrays. - bSector(); +namespace magneticfield { - /// Constructor from list of volumes - bSector(handles::const_iterator begin, handles::const_iterator end); + class bSector { + public: + /// Default ctor is needed to have arrays. + bSector(); - /// Destructor - ~bSector(); + /// Constructor from list of volumes + bSector(handles::const_iterator begin, handles::const_iterator end, bool debugVal = false); - /// Distance from center along normal of sectors. - const float RN() const { return volumes.front()->RN(); } + /// Destructor + ~bSector() = default; - /// Return all volumes in this sector - const handles& getVolumes() const { return volumes; } + /// Distance from center along normal of sectors. + const float RN() const { return volumes.front()->RN(); } - /// Construct the MagBSector upon request. - MagBSector* buildMagBSector() const; + /// Return all volumes in this sector + const handles& getVolumes() const { return volumes; } + + /// Construct the MagBSector upon request. + MagBSector* buildMagBSector() const; + + private: + std::vector rods; // the rods in this layer + handles volumes; // pointers to all volumes in the sector + mutable MagBSector* msector; + bool debug; // Allow assignment from other bSector objects + }; +} // namespace magneticfield -private: - std::vector rods; // the rods in this layer - handles volumes; // pointers to all volumes in the sector - mutable MagBSector* msector; -}; #endif diff --git a/MagneticField/GeomBuilder/src/bSlab.cc b/MagneticField/GeomBuilder/src/bSlab.cc index fd8b9b109a8be..2e526b3785fe9 100644 --- a/MagneticField/GeomBuilder/src/bSlab.cc +++ b/MagneticField/GeomBuilder/src/bSlab.cc @@ -13,16 +13,15 @@ using namespace SurfaceOrientation; using namespace std; +using namespace magneticfield; -MagGeoBuilderFromDDD::bSlab::~bSlab() {} - -MagGeoBuilderFromDDD::bSlab::bSlab(handles::const_iterator begin, handles::const_iterator end) - : volumes(begin, end), mslab(nullptr) { +bSlab::bSlab(handles::const_iterator begin, handles::const_iterator end, bool debugVal) + : volumes(begin, end), mslab(nullptr), debug(debugVal) { if (volumes.size() > 1) { // Sort volumes by dphi i.e. phi(j)-phi(i) > 0 if j>1. precomputed_value_sort(volumes.begin(), volumes.end(), ExtractPhiMax(), LessDPhi()); - if (MagGeoBuilderFromDDD::debug) + if (debug) cout << " Slab has " << volumes.size() << " volumes" << endl; // Check that all volumes have the same dZ @@ -33,7 +32,7 @@ MagGeoBuilderFromDDD::bSlab::bSlab(handles::const_iterator begin, handles::const const float epsilon = 0.001; if (fabs(Zmax - (*i)->surface(zplus).position().z()) > epsilon || fabs(Zmin - (*i)->surface(zminus).position().z()) > epsilon) { - if (MagGeoBuilderFromDDD::debug) + if (debug) cout << "*** WARNING: slabs Z coords not matching: D_Zmax = " << fabs(Zmax - (*i)->surface(zplus).position().z()) << " D_Zmin = " << fabs(Zmin - (*i)->surface(zminus).position().z()) << endl; @@ -42,11 +41,11 @@ MagGeoBuilderFromDDD::bSlab::bSlab(handles::const_iterator begin, handles::const } } -Geom::Phi MagGeoBuilderFromDDD::bSlab::minPhi() const { return volumes.front()->minPhi(); } +Geom::Phi bSlab::minPhi() const { return volumes.front()->minPhi(); } -Geom::Phi MagGeoBuilderFromDDD::bSlab::maxPhi() const { return volumes.back()->maxPhi(); } +Geom::Phi bSlab::maxPhi() const { return volumes.back()->maxPhi(); } -MagBSlab* MagGeoBuilderFromDDD::bSlab::buildMagBSlab() const { +MagBSlab* bSlab::buildMagBSlab() const { if (mslab == nullptr) { vector mVols; for (handles::const_iterator vol = volumes.begin(); vol != volumes.end(); ++vol) { diff --git a/MagneticField/GeomBuilder/src/bSlab.h b/MagneticField/GeomBuilder/src/bSlab.h index 73f539b4d2625..69969d3ac1c1a 100644 --- a/MagneticField/GeomBuilder/src/bSlab.h +++ b/MagneticField/GeomBuilder/src/bSlab.h @@ -9,37 +9,51 @@ * \author N. Amapane - INFN Torino */ -#include "MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.h" -#include "MagneticField/GeomBuilder/src/volumeHandle.h" -#include "DataFormats/GeometryVector/interface/Pi.h" +#include "DataFormats/GeometryVector/interface/Phi.h" + +#include "MagneticField/VolumeGeometry/interface/MagVolume6Faces.h" +#include "MagneticField/Layers/interface/MagBSlab.h" +#include "MagneticField/Layers/interface/MagVerbosity.h" +#include "MagneticField/GeomBuilder/src/BaseVolumeHandle.h" + +#include "Utilities/General/interface/precomputed_value_sort.h" + +#include +#include class MagBSlab; -class MagGeoBuilderFromDDD::bSlab { -public: - /// Constructor from list of volumes - bSlab(handles::const_iterator begin, handles::const_iterator end); +namespace magneticfield { + + using std::vector; + + class bSlab { + public: + /// Constructor from list of volumes + bSlab(handles::const_iterator startIter, handles::const_iterator endIter, bool debugVal = false); - /// Destructor - ~bSlab(); + /// Destructor + ~bSlab() = default; - /// Distance from center along sector normal. - const float RN() const { return volumes.front()->RN(); } + /// Distance from center along sector normal. + const float RN() const { return volumes.front()->RN(); } - /// Boundary in phi. - // FIXME: use volumeHandle [max|min]Phi, which returns phi at median of - // phi plane (not absolute limits). Used by: bRod ctor (only for dphi) - Geom::Phi minPhi() const; + /// Boundary in phi. + // FIXME: use volumeHandle [max|min]Phi, which returns phi at median of + // phi plane (not absolute limits). Used by: bRod ctor (only for dphi) + Geom::Phi minPhi() const; - /// Boundary in phi. - Geom::Phi maxPhi() const; + /// Boundary in phi. + Geom::Phi maxPhi() const; - /// Construct the MagBSlab upon request. - MagBSlab* buildMagBSlab() const; + /// Construct the MagBSlab upon request. + MagBSlab* buildMagBSlab() const; -private: - handles volumes; // pointers to all volumes in the slab - mutable MagBSlab* mslab; -}; + private: + handles volumes; // pointers to all volumes in the slab + mutable MagBSlab* mslab; + bool debug; // Allow assignment + }; +} // namespace magneticfield #endif diff --git a/MagneticField/GeomBuilder/src/buildBox.icc b/MagneticField/GeomBuilder/src/buildBox.icc index 9b7a0ea768c95..53b372511ef5e 100644 --- a/MagneticField/GeomBuilder/src/buildBox.icc +++ b/MagneticField/GeomBuilder/src/buildBox.icc @@ -4,14 +4,15 @@ * \author N. Amapane - INFN Torino */ -void MagGeoBuilderFromDDD::volumeHandle::buildBox(const DDExpandedView &fv) { - if (MagGeoBuilderFromDDD::debug) - cout << "Building box surfaces...: " << endl; +void volumeHandle::buildBox() { + LogDebug("magneticfield::volumeHandle") << "Building box surfaces...: "; DDBox box(solid); - double halfX = box.halfX() / cm; - double halfY = box.halfY() / cm; - double halfZ = box.halfZ() / cm; + // Old DD needs mm to cm conversion, but DD4hep needs no conversion. + // convertUnits should be defined appropriately. + double halfX = convertUnits(box.halfX()); + double halfY = convertUnits(box.halfY()); + double halfZ = convertUnits(box.halfZ()); // Global vectors of the normals to X, Y, Z axes GlobalVector planeXAxis = refPlane->toGlobal(LocalVector(1, 0, 0)); @@ -39,9 +40,9 @@ void MagGeoBuilderFromDDD::volumeHandle::buildBox(const DDExpandedView &fv) { Surface::RotationType rot_phi; Surface::RotationType rot_Z = Surface::RotationType(planeXAxis, planeYAxis); - if (fabs(rnX) > fabs(rnY)) { + if (std::abs(rnX) > std::abs(rnY)) { // X is ~parallel to global R dir, Y is along +/- phi - theRN = fabs(rnX); + theRN = std::abs(rnX); if (rnX < 0) { halfX = -halfX; halfY = -halfY; @@ -55,7 +56,7 @@ void MagGeoBuilderFromDDD::volumeHandle::buildBox(const DDExpandedView &fv) { rot_phi = Surface::RotationType(planeZAxis, planeXAxis); // opposite to y axis } else { // Y is ~parallel to global R dir, X is along +/- phi - theRN = fabs(rnY); + theRN = std::abs(rnY); if (rnY < 0) { halfX = -halfX; halfY = -halfY; @@ -69,28 +70,27 @@ void MagGeoBuilderFromDDD::volumeHandle::buildBox(const DDExpandedView &fv) { rot_phi = Surface::RotationType(planeZAxis, planeYAxis); // opposite to x axis } - if (MagGeoBuilderFromDDD::debug) - cout << " halfX: " << halfX << " halfY: " << halfY << " halfZ: " << halfZ << endl - << "RN " << theRN << endl; + LogDebug("magneticfield::volumeHandle") << " halfX: " << halfX << " halfY: " << halfY << " halfZ: " << halfZ << endl + << "RN " << theRN << endl; - if (MagGeoBuilderFromDDD::debug) - cout << "pos_outer " << pos_outer << " " << pos_outer.perp() << " " << pos_outer.phi() << endl - << "pos_inner " << pos_inner << " " << pos_inner.perp() << " " << pos_inner.phi() << endl - << "pos_zplus " << pos_zplus << " " << pos_zplus.perp() << " " << pos_zplus.phi() << endl - << "pos_zminus " << pos_zminus << " " << pos_zminus.perp() << " " << pos_zminus.phi() << endl - << "pos_phiplus " << pos_phiplus << " " << pos_phiplus.perp() << " " << pos_phiplus.phi() << endl - << "pos_phiminus " << pos_phiminus << " " << pos_phiminus.perp() << " " << pos_phiminus.phi() << endl; + LogDebug("magneticfield::volumeHandle") + << "pos_outer " << pos_outer << " " << pos_outer.perp() << " " << pos_outer.phi() << endl + << "pos_inner " << pos_inner << " " << pos_inner.perp() << " " << pos_inner.phi() << endl + << "pos_zplus " << pos_zplus << " " << pos_zplus.perp() << " " << pos_zplus.phi() << endl + << "pos_zminus " << pos_zminus << " " << pos_zminus.perp() << " " << pos_zminus.phi() << endl + << "pos_phiplus " << pos_phiplus << " " << pos_phiplus.perp() << " " << pos_phiplus.phi() << endl + << "pos_phiminus " << pos_phiminus << " " << pos_phiminus.perp() << " " << pos_phiminus.phi() << endl; // Check ordering. - if (MagGeoBuilderFromDDD::debug) { + if (debug) { if (pos_outer.perp() < pos_inner.perp()) { - cout << "*** WARNING: pos_outer < pos_inner for box" << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pos_outer < pos_inner for box" << endl; } if (pos_zplus.z() < pos_zminus.z()) { - cout << "*** WARNING: pos_zplus < pos_zminus for box" << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pos_zplus < pos_zminus for box" << endl; } if (Geom::Phi(pos_phiplus.phi() - pos_phiminus.phi()) < 0.) { - cout << "*** WARNING: pos_phiplus < pos_phiminus for box" << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pos_phiplus < pos_phiminus for box" << endl; } } @@ -102,15 +102,13 @@ void MagGeoBuilderFromDDD::volumeHandle::buildBox(const DDExpandedView &fv) { surfaces[phiplus] = new Plane(pos_phiplus, rot_phi); surfaces[phiminus] = new Plane(pos_phiminus, rot_phi); - if (MagGeoBuilderFromDDD::debug) { - cout << "rot_R " << surfaces[outer]->toGlobal(LocalVector(0., 0., 1.)) << endl - << "rot_Z " << surfaces[zplus]->toGlobal(LocalVector(0., 0., 1.)) << endl - << "rot_phi " << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)) << endl; - } + LogDebug("magneticfield::volumeHandle") << "rot_R " << surfaces[outer]->toGlobal(LocalVector(0., 0., 1.)) << endl + << "rot_Z " << surfaces[zplus]->toGlobal(LocalVector(0., 0., 1.)) << endl + << "rot_phi " << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)); // Save volume boundaries - theRMin = fabs(surfaces[inner]->toLocal(GlobalPoint(0, 0, 0)).z()); - theRMax = fabs(surfaces[outer]->toLocal(GlobalPoint(0, 0, 0)).z()); + theRMin = std::abs(surfaces[inner]->toLocal(GlobalPoint(0, 0, 0)).z()); + theRMax = std::abs(surfaces[outer]->toLocal(GlobalPoint(0, 0, 0)).z()); // FIXME: use phi of middle plane of phiminus surface. Is not the absolute phimin! thePhiMin = surfaces[phiminus]->position().phi(); } diff --git a/MagneticField/GeomBuilder/src/buildCons.icc b/MagneticField/GeomBuilder/src/buildCons.icc index 23f70af2286cd..3e38a809ee94b 100644 --- a/MagneticField/GeomBuilder/src/buildCons.icc +++ b/MagneticField/GeomBuilder/src/buildCons.icc @@ -6,28 +6,28 @@ #include "DataFormats/GeometrySurface/interface/SimpleConeBounds.h" -void MagGeoBuilderFromDDD::volumeHandle::buildCons(const DDExpandedView &fv) { - if (MagGeoBuilderFromDDD::debug) - cout << "Building cons surfaces...: " << endl; +void volumeHandle::buildCons() { + LogDebug("magneticfield::volumeHandle") << "Building cons surfaces...: " << endl; DDCons cons(solid); - double zhalf = cons.zhalf() / cm; - double rInMinusZ = cons.rInMinusZ() / cm; - double rOutMinusZ = cons.rOutMinusZ() / cm; - double rInPlusZ = cons.rInPlusZ() / cm; - double rOutPlusZ = cons.rOutPlusZ() / cm; + // Old DD needs mm to cm conversion, but DD4hep needs no conversion. + // convertUnits should be defined appropriately. + double zhalf = convertUnits(cons.zhalf()); + double rInMinusZ = convertUnits(cons.rInMinusZ()); + double rOutMinusZ = convertUnits(cons.rOutMinusZ()); + double rInPlusZ = convertUnits(cons.rInPlusZ()); + double rOutPlusZ = convertUnits(cons.rOutPlusZ()); double startPhi = cons.phiFrom(); double deltaPhi = cons.deltaPhi(); - if (MagGeoBuilderFromDDD::debug) - cout << "zhalf " << zhalf << endl - << "rInMinusZ " << rInMinusZ << endl - << "rOutMinusZ " << rOutMinusZ << endl - << "rInPlusZ " << rInPlusZ << endl - << "rOutPlusZ " << rOutPlusZ << endl - << "phiFrom " << startPhi << endl - << "deltaPhi " << deltaPhi << endl; + LogDebug("magneticfield::volumeHandle") << "zhalf " << zhalf << endl + << "rInMinusZ " << rInMinusZ << endl + << "rOutMinusZ " << rOutMinusZ << endl + << "rInPlusZ " << rInPlusZ << endl + << "rOutPlusZ " << rOutPlusZ << endl + << "phiFrom " << startPhi << endl + << "deltaPhi " << deltaPhi << endl; // recalculate center: (for a DDCons, DDD gives 0,0,Z) double rZmin = (rInMinusZ + rOutMinusZ) / 2.; @@ -40,7 +40,7 @@ void MagGeoBuilderFromDDD::volumeHandle::buildCons(const DDExpandedView &fv) { const double epsilon = 1e-5; - if (abs(rInPlusZ - rInMinusZ) < epsilon) { // Cylinder + if (std::abs(rInPlusZ - rInMinusZ) < epsilon) { // Cylinder // FIXME: use builder surfaces[inner] = new Cylinder(rInMinusZ, Surface::PositionType(), Surface::RotationType()); @@ -51,10 +51,8 @@ void MagGeoBuilderFromDDD::volumeHandle::buildCons(const DDExpandedView &fv) { surfaces[inner] = new Cone(Surface::PositionType(0, 0, center_.z()), Surface::RotationType(), cb.vertex(), cb.openingAngle()); } - - if (abs(rOutPlusZ - rOutMinusZ) < epsilon) { // Cylinder + if (std::abs(rOutPlusZ - rOutMinusZ) < epsilon) { // Cylinder surfaces[outer] = new Cylinder(rOutMinusZ, Surface::PositionType(0, 0, center_.z()), Surface::RotationType()); - } else { // Cone // FIXME: trick to compute vertex and angle... SimpleConeBounds cb(center_.z() - zhalf, rOutMinusZ, rOutMinusZ, center_.z() + zhalf, rOutPlusZ, rOutPlusZ); @@ -62,22 +60,12 @@ void MagGeoBuilderFromDDD::volumeHandle::buildCons(const DDExpandedView &fv) { surfaces[outer] = new Cone(Surface::PositionType(0, 0, center_.z()), Surface::RotationType(), cb.vertex(), cb.openingAngle()); - if (MagGeoBuilderFromDDD::debug) - cout << "Outer surface: cone, vtx: " << cb.vertex() << " angle " << cb.openingAngle() << endl; - - // cout << (cb.vertex().z()-center_.z()-zhalf) *tan(cb.openingAngle()) << endl; - // cout << (cb.vertex().z()-center_.z()+zhalf) *tan(cb.openingAngle()) << endl; + LogDebug("magneticfield::volumeHandle") + << "Outer surface: cone, vtx: " << cb.vertex() << " angle " << cb.openingAngle() << endl; } - // All other surfaces buildPhiZSurf(startPhi, deltaPhi, zhalf, rCentr); - // Check ordering. //FIXME - // if (dynamic_cast(&(*surfaces[outer]))->bounds().width() < - // dynamic_cast(&(*surfaces[inner]))->bounds().width()){ - // cout << "*** WARNING: pos_outer < pos_inner " << endl; - // } - // Save volume boundaries theRMin = min(rInMinusZ, rInPlusZ); theRMax = max(rOutMinusZ, rOutPlusZ); diff --git a/MagneticField/GeomBuilder/src/buildPseudoTrap.icc b/MagneticField/GeomBuilder/src/buildPseudoTrap.icc index 3f3ff10376290..6e8978e8521f3 100644 --- a/MagneticField/GeomBuilder/src/buildPseudoTrap.icc +++ b/MagneticField/GeomBuilder/src/buildPseudoTrap.icc @@ -1,5 +1,6 @@ /* * Compute parameters for a pseudotrapezoid . + * NOTE: The pseudotrapezoid is no longer needed or support for DD4hep. * For pseudotraps the refPlane is parallel to beam line; transformation: * * Global(for vol at pi/2) Local @@ -10,42 +11,41 @@ * \author N. Amapane - INFN Torino */ -void MagGeoBuilderFromDDD::volumeHandle::buildPseudoTrap(const DDExpandedView &fv) { - if (MagGeoBuilderFromDDD::debug) - cout << "Building PseudoTrap surfaces...: " << endl; +void volumeHandle::buildPseudoTrap() { + LogDebug("magneticfield::volumeHandle") << "Building PseudoTrap surfaces...: " << endl; DDPseudoTrap ptrap(solid); - double halfZ = ptrap.halfZ() / cm; - double x1 = ptrap.x1() / cm; - double x2 = ptrap.x2() / cm; - double y1 = ptrap.y1() / cm; - double y2 = ptrap.y2() / cm; - double radius = ptrap.radius() / cm; + // Old DD needs mm to cm conversion, but DD4hep needs no conversion. + // convertUnits should be defined appropriately. + double halfZ = convertUnits(ptrap.halfZ()); + double x1 = convertUnits(ptrap.x1()); + double x2 = convertUnits(ptrap.x2()); + double y1 = convertUnits(ptrap.y1()); + double y2 = convertUnits(ptrap.y2()); + double radius = convertUnits(ptrap.radius()); bool atMinusZ = ptrap.atMinusZ(); - if (MagGeoBuilderFromDDD::debug) - cout << "halfZ " << halfZ << endl - << "x1 " << x1 << endl - << "x2 " << x2 << endl - << "y1 " << y1 << endl - << "y2 " << y2 << endl - << "radius " << radius << endl - << "atMinusZ " << atMinusZ << endl; + LogDebug("magneticfield::volumeHandle") << "halfZ " << halfZ << endl + << "x1 " << x1 << endl + << "x2 " << x2 << endl + << "y1 " << y1 << endl + << "y2 " << y2 << endl + << "radius " << radius << endl + << "atMinusZ " << atMinusZ << endl; // Just check assumptions on parameters... const double epsilon = 1e-5; - if (MagGeoBuilderFromDDD::debug) { + if (debug) { if (y1 - y2 > epsilon) { - cout << "*** WARNING: unexpected pseudotrapezoid parameters." << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: unexpected pseudotrapezoid parameters." << endl; } // Check that volume is convex (no concave volume in current geometry...) if (radius * (atMinusZ ? -1. : +1) < 0.) { - cout << "*** WARNING: pseudotrapezoid is concave" << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pseudotrapezoid is concave" << endl; } } - // CAVEAT: pseudotraps are rotated in a different way as traps, // since they have local Z along global R... GlobalVector planeXAxis = refPlane->toGlobal(LocalVector(1, 0, 0)); @@ -55,9 +55,9 @@ void MagGeoBuilderFromDDD::volumeHandle::buildPseudoTrap(const DDExpandedView &f //FIXME assumption: here we do use the assumption on the orientation / //of local Z (see above) GlobalVector Rvol(refPlane->position().x(), refPlane->position().y(), refPlane->position().z()); - theRN = fabs(planeZAxis.dot(Rvol)); + theRN = std::abs(planeZAxis.dot(Rvol)); - double fR = fabs(radius); + double fR = std::abs(radius); Sides cyl_side; Sides plane_side; if (atMinusZ) { @@ -81,42 +81,41 @@ void MagGeoBuilderFromDDD::volumeHandle::buildPseudoTrap(const DDExpandedView &f else rcheck = refPlane->toGlobal(LocalPoint(x2, 0., +halfZ)).perp(); - if (MagGeoBuilderFromDDD::debug) { - if (fabs(rcheck - fR) > 100. * epsilon) { //FIXME! - cout << setprecision(10); - cout << "*** WARNING: Cylinder surface not centered on beam axis " << rcheck << " " << fR << " " - << fabs(rcheck - fR) << endl; - cout << setprecision(6); + if (debug) { + if (std::abs(rcheck - fR) > 100. * epsilon) { //FIXME! + LogDebug("magneticfield::volumeHandle") << setprecision(10); + LogDebug("magneticfield::volumeHandle") << "*** WARNING: Cylinder surface not centered on beam axis " << rcheck + << " " << fR << " " << std::abs(rcheck - fR) << endl; + LogDebug("magneticfield::volumeHandle") << setprecision(6); } } - if (MagGeoBuilderFromDDD::debug) { - cout << "RN " << theRN << endl - << "pos_Rplane " << pos_Rplane << " " << pos_Rplane.perp() << " " << pos_Rplane.phi() << endl - << "pos_zplus " << pos_zplus << " " << pos_zplus.perp() << " " << pos_zplus.phi() << endl - << "pos_zminus " << pos_zminus << " " << pos_zminus.perp() << " " << pos_zminus.phi() << endl - << "pos_phiplus " << pos_phiplus << " " << pos_phiplus.perp() << " " << pos_phiplus.phi() << endl - << "pos_phiminus " << pos_phiminus << " " << pos_phiminus.perp() << " " << pos_phiminus.phi() << endl; + if (debug) { + LogDebug("magneticfield::volumeHandle") + << "RN " << theRN << endl + << "pos_Rplane " << pos_Rplane << " " << pos_Rplane.perp() << " " << pos_Rplane.phi() << endl + << "pos_zplus " << pos_zplus << " " << pos_zplus.perp() << " " << pos_zplus.phi() << endl + << "pos_zminus " << pos_zminus << " " << pos_zminus.perp() << " " << pos_zminus.phi() << endl + << "pos_phiplus " << pos_phiplus << " " << pos_phiplus.perp() << " " << pos_phiplus.phi() << endl + << "pos_phiminus " << pos_phiminus << " " << pos_phiminus.perp() << " " << pos_phiminus.phi() << endl; } // Check ordering. if ((pos_Rplane.perp() < radius) == atMinusZ) { - cout << "*** WARNING: pos_outer < pos_inner for pseudotrapezoid" << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pos_outer < pos_inner for pseudotrapezoid" << endl; } if (pos_zplus.z() < pos_zminus.z()) { - cout << "*** WARNING: pos_zplus < pos_zminus for pseudotrapezoid" << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pos_zplus < pos_zminus for pseudotrapezoid" << endl; } if (Geom::Phi(pos_phiplus.phi() - pos_phiminus.phi()) < 0.) { - cout << "*** WARNING: pos_phiplus < pos_phiminus for pseudotrapezoid" << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pos_phiplus < pos_phiminus for pseudotrapezoid" << endl; } GlobalVector z_phiplus = (refPlane->toGlobal(LocalVector((x2 - x1) / 2., 0., halfZ))).unit(); GlobalVector z_phiminus = (refPlane->toGlobal(LocalVector(-(x2 - x1) / 2., 0., halfZ))).unit(); - if (MagGeoBuilderFromDDD::debug) { - cout << "z_phiplus " << z_phiplus << " " << z_phiplus.phi() << endl; - cout << "z_phiminus " << z_phiminus << " " << z_phiminus.phi() << endl; - } + LogDebug("magneticfield::volumeHandle") << "z_phiplus " << z_phiplus << " " << z_phiplus.phi() << endl; + LogDebug("magneticfield::volumeHandle") << "z_phiminus " << z_phiminus << " " << z_phiminus.phi() << endl; Surface::RotationType rot_R(planeYAxis, planeXAxis); Surface::RotationType rot_Z(planeXAxis, planeZAxis); @@ -132,17 +131,16 @@ void MagGeoBuilderFromDDD::volumeHandle::buildPseudoTrap(const DDExpandedView &f surfaces[phiplus] = new Plane(pos_phiplus, rot_phiplus); surfaces[phiminus] = new Plane(pos_phiminus, rot_phiminus); - if (MagGeoBuilderFromDDD::debug) { - cout << "rot_R " << surfaces[plane_side]->toGlobal(LocalVector(0., 0., 1.)) << endl - << "rot_Z " << surfaces[zplus]->toGlobal(LocalVector(0., 0., 1.)) << endl - << "rot_phi+ " << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " - << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl - << "rot_phi- " << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " - << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl; - } + LogDebug("magneticfield::volumeHandle") + << "rot_R " << surfaces[plane_side]->toGlobal(LocalVector(0., 0., 1.)) << endl + << "rot_Z " << surfaces[zplus]->toGlobal(LocalVector(0., 0., 1.)) << endl + << "rot_phi+ " << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " + << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl + << "rot_phi- " << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " + << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)).phi(); // Save volume boundaries - double R1 = fabs(surfaces[plane_side]->toLocal(GlobalPoint(0, 0, 0)).z()); + double R1 = std::abs(surfaces[plane_side]->toLocal(GlobalPoint(0, 0, 0)).z()); theRMin = min(fR, R1); theRMax = max(fR, R1); // FIXME: use phi of middle plane of phiminus surface. Is not the absolute phimin! diff --git a/MagneticField/GeomBuilder/src/buildTrap.icc b/MagneticField/GeomBuilder/src/buildTrap.icc index 8f2e728676cfe..4935da3e691c6 100644 --- a/MagneticField/GeomBuilder/src/buildTrap.icc +++ b/MagneticField/GeomBuilder/src/buildTrap.icc @@ -4,43 +4,43 @@ * \author N. Amapane - INFN Torino */ -void MagGeoBuilderFromDDD::volumeHandle::buildTrap(const DDExpandedView &fv) { - if (MagGeoBuilderFromDDD::debug) - cout << "Building trapezoid surfaces...: " << endl; +void volumeHandle::buildTrap() { + LogDebug("magneticfield::volumeHandle") << "Building trapezoid surfaces...: " << endl; DDTrap trap(solid); //FIXME - double x1 = trap.x1() / cm; - double x2 = trap.x2() / cm; - double x3 = trap.x3() / cm; - double x4 = trap.x4() / cm; - double y1 = trap.y1() / cm; - double y2 = trap.y2() / cm; + // Old DD needs mm to cm conversion, but DD4hep needs no conversion. + // convertUnits should be defined appropriately. + double x1 = convertUnits(trap.x1()); + double x2 = convertUnits(trap.x2()); + double x3 = convertUnits(trap.x3()); + double x4 = convertUnits(trap.x4()); + double y1 = convertUnits(trap.y1()); + double y2 = convertUnits(trap.y2()); double theta = trap.theta(); double phi = trap.phi(); - double halfZ = trap.halfZ() / cm; + double halfZ = convertUnits(trap.halfZ()); double alpha1 = trap.alpha1(); double alpha2 = trap.alpha2(); - if (MagGeoBuilderFromDDD::debug) - cout << "x1 " << x1 << endl - << "x2 " << x2 << endl - << "x3 " << x3 << endl - << "x4 " << x4 << endl - << "y1 " << y1 << endl - << "y2 " << y2 << endl - << "theta " << theta << endl - << "phi " << phi << endl - << "halfZ " << halfZ << endl - << "alpha1 " << alpha1 << endl - << "alpha2 " << alpha2 << endl; + LogDebug("magneticfield::volumeHandle") << "x1 " << x1 << endl + << "x2 " << x2 << endl + << "x3 " << x3 << endl + << "x4 " << x4 << endl + << "y1 " << y1 << endl + << "y2 " << y2 << endl + << "theta " << theta << endl + << "phi " << phi << endl + << "halfZ " << halfZ << endl + << "alpha1 " << alpha1 << endl + << "alpha2 " << alpha2 << endl; // Just check assumptions on parameters... - if (MagGeoBuilderFromDDD::debug) { + if (debug) { const double epsilon = 1e-5; if (theta > epsilon || phi > epsilon || y1 - y2 > epsilon || x1 - x3 > epsilon || x2 - x4 > epsilon || alpha1 - alpha2 > epsilon) { - cout << "*** WARNING: unexpected trapezoid parameters." << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: unexpected trapezoid parameters." << endl; } } @@ -59,7 +59,7 @@ void MagGeoBuilderFromDDD::volumeHandle::buildTrap(const DDExpandedView &fv) { // FIXME Assumption: it is assumed that local Y is always along global R // (true for version 1103l, not necessarily in the future) GlobalVector Rvol(refPlane->position().x(), refPlane->position().y(), refPlane->position().z()); - theRN = fabs(planeYAxis.dot(Rvol)); + theRN = std::abs(planeYAxis.dot(Rvol)); if (planeZAxis.z() < 0.) { x1 = -x1; @@ -76,25 +76,25 @@ void MagGeoBuilderFromDDD::volumeHandle::buildTrap(const DDExpandedView &fv) { GlobalPoint pos_phiplus(refPlane->toGlobal(LocalPoint(-halfX, 0., 0.))); GlobalPoint pos_phiminus(refPlane->toGlobal(LocalPoint(halfX, 0., 0.))); - if (MagGeoBuilderFromDDD::debug) - cout << "RN " << theRN << endl - << "pos_outer " << pos_outer << " " << pos_outer.perp() << " " << pos_outer.phi() << endl - << "pos_inner " << pos_inner << " " << pos_inner.perp() << " " << pos_inner.phi() << endl - << "pos_zplus " << pos_zplus << " " << pos_zplus.perp() << " " << pos_zplus.phi() << endl - << "pos_zminus " << pos_zminus << " " << pos_zminus.perp() << " " << pos_zminus.phi() << endl - << "pos_phiplus " << pos_phiplus << " " << pos_phiplus.perp() << " " << pos_phiplus.phi() << endl - << "pos_phiminus " << pos_phiminus << " " << pos_phiminus.perp() << " " << pos_phiminus.phi() << endl; + LogDebug("magneticfield::volumeHandle") + << "RN " << theRN << endl + << "pos_outer " << pos_outer << " " << pos_outer.perp() << " " << pos_outer.phi() << endl + << "pos_inner " << pos_inner << " " << pos_inner.perp() << " " << pos_inner.phi() << endl + << "pos_zplus " << pos_zplus << " " << pos_zplus.perp() << " " << pos_zplus.phi() << endl + << "pos_zminus " << pos_zminus << " " << pos_zminus.perp() << " " << pos_zminus.phi() << endl + << "pos_phiplus " << pos_phiplus << " " << pos_phiplus.perp() << " " << pos_phiplus.phi() << endl + << "pos_phiminus " << pos_phiminus << " " << pos_phiminus.perp() << " " << pos_phiminus.phi() << endl; // Check ordering. - if (MagGeoBuilderFromDDD::debug) { + if (debug) { if (pos_outer.perp() < pos_inner.perp()) { - cout << "*** WARNING: pos_outer < pos_inner for trapezoid" << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pos_outer < pos_inner for trapezoid" << endl; } if (pos_zplus.z() < pos_zminus.z()) { - cout << "*** WARNING: pos_zplus < pos_zminus for trapezoid" << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pos_zplus < pos_zminus for trapezoid" << endl; } if (Geom::Phi(pos_phiplus.phi() - pos_phiminus.phi()) < 0.) { - cout << "*** WARNING: pos_phiplus < pos_phiminus for trapezoid" << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pos_phiplus < pos_phiminus for trapezoid" << endl; } } @@ -103,10 +103,8 @@ void MagGeoBuilderFromDDD::volumeHandle::buildTrap(const DDExpandedView &fv) { GlobalVector y_phiminusV = (refPlane->toGlobal(LocalVector((tan(alpha1) * y1 + (x2 - x1) / 2.), y1, 0.))); GlobalVector y_phiminus = y_phiminusV.unit(); - if (MagGeoBuilderFromDDD::debug) { - cout << "y_phiplus " << y_phiplus << endl; - cout << "y_phiminus " << y_phiminus << endl; - } + LogDebug("magneticfield::volumeHandle") << "y_phiplus " << y_phiplus << endl; + LogDebug("magneticfield::volumeHandle") << "y_phiminus " << y_phiminus << endl; Surface::RotationType rot_R(planeZAxis, planeXAxis); Surface::RotationType rot_Z(planeXAxis, planeYAxis); @@ -122,8 +120,8 @@ void MagGeoBuilderFromDDD::volumeHandle::buildTrap(const DDExpandedView &fv) { surfaces[phiminus] = new Plane(pos_phiminus, rot_phiminus); // Save volume boundaries - theRMin = fabs(surfaces[inner]->toLocal(GlobalPoint(0, 0, 0)).z()); - theRMax = fabs(surfaces[outer]->toLocal(GlobalPoint(0, 0, 0)).z()); + theRMin = std::abs(surfaces[inner]->toLocal(GlobalPoint(0, 0, 0)).z()); + theRMax = std::abs(surfaces[outer]->toLocal(GlobalPoint(0, 0, 0)).z()); // Setting "absolute" phi boundary spots a problem in V85l // e.g. vol 139 (origin: small inconsistency in volumes.txt) @@ -135,16 +133,17 @@ void MagGeoBuilderFromDDD::volumeHandle::buildTrap(const DDExpandedView &fv) { thePhiMin = min((pos_phiminus + y_phiminusV).phi(), (pos_phiminus - y_phiminusV).phi()); } - if (MagGeoBuilderFromDDD::debug) { - cout << "rot_R " << surfaces[outer]->toGlobal(LocalVector(0., 0., 1.)) << endl - << "rot_Z " << surfaces[zplus]->toGlobal(LocalVector(0., 0., 1.)) << endl - << "rot_phi+ " << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " - << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl - << "rot_phi- " << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " - << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl; - - if (fabs(surfaces[phiminus]->position().phi() - thePhiMin) > 1e-5) { - cout << " Trapez. phiMin = " << thePhiMin << endl; + if (debug) { + LogDebug("magneticfield::volumeHandle") + << "rot_R " << surfaces[outer]->toGlobal(LocalVector(0., 0., 1.)) << endl + << "rot_Z " << surfaces[zplus]->toGlobal(LocalVector(0., 0., 1.)) << endl + << "rot_phi+ " << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " + << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl + << "rot_phi- " << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " + << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl; + + if (std::abs(surfaces[phiminus]->position().phi() - thePhiMin) > 1e-5) { + LogDebug("magneticfield::volumeHandle") << " Trapez. phiMin = " << thePhiMin << endl; } } } diff --git a/MagneticField/GeomBuilder/src/buildTruncTubs.icc b/MagneticField/GeomBuilder/src/buildTruncTubs.icc index 61603cc305fe9..a19ec4d74d399 100644 --- a/MagneticField/GeomBuilder/src/buildTruncTubs.icc +++ b/MagneticField/GeomBuilder/src/buildTruncTubs.icc @@ -9,30 +9,26 @@ * \author N. Amapane - INFN Torino */ -void MagGeoBuilderFromDDD::volumeHandle::buildTruncTubs(const DDExpandedView &fv) { - if (MagGeoBuilderFromDDD::debug) - cout << "Building TruncTubs surfaces...: " << endl; +void volumeHandle::buildTruncTubs() { + LogDebug("magneticfield::volumeHandle") << "Building TruncTubs surfaces...: " << endl; DDTruncTubs tubs(solid); - double zhalf = tubs.zHalf() / cm; // half of the z-Axis - double rIn = tubs.rIn() / cm; // inner radius - double rOut = tubs.rOut() / cm; // outer radius - double startPhi = tubs.startPhi(); // angular start of the tube-section - double deltaPhi = tubs.deltaPhi(); // angular span of the tube-section - double cutAtStart = tubs.cutAtStart() / cm; // truncation at begin of the tube-section - double cutAtDelta = tubs.cutAtDelta() / cm; // truncation at end of the tube-section - bool cutInside = tubs.cutInside(); // true, if truncation is on the inner side of the tube-section + // Old DD needs mm to cm conversion, but DD4hep needs no conversion. + // convertUnits should be defined appropriately. + double zhalf = convertUnits(tubs.zHalf()); // half of the z-Axis + double rIn = convertUnits(tubs.rIn()); // inner radius + double rOut = convertUnits(tubs.rOut()); // outer radius + double startPhi = tubs.startPhi(); // angular start of the tube-section + double deltaPhi = tubs.deltaPhi(); // angular span of the tube-section + double cutAtStart = convertUnits(tubs.cutAtStart()); // truncation at begin of the tube-section + double cutAtDelta = convertUnits(tubs.cutAtDelta()); // truncation at end of the tube-section + bool cutInside = tubs.cutInside(); // true, if truncation is on the inner side of the tube-section - if (MagGeoBuilderFromDDD::debug) - cout << "zhalf " << zhalf << endl - << "rIn " << rIn << endl - << "rOut " << rOut << endl - << "startPhi " << startPhi << endl - << "deltaPhi " << deltaPhi << endl - << "cutAtStart " << cutAtStart << endl - << "cutAtDelta " << cutAtDelta << endl - << "cutInside " << cutInside << endl; + LogDebug("magneticfield::volumeHandle") + << " zhalf " << zhalf << "\n rIn " << rIn << "\n rOut " << rOut << "\n startPhi " << startPhi + << "\n deltaPhi " << deltaPhi << "\n cutAtStart " << cutAtStart << "\n cutAtDelta " << cutAtDelta + << "\n cutInside " << cutInside; // // GlobalVector planeXAxis = refPlane->toGlobal(LocalVector( 1, 0, 0)); @@ -89,18 +85,18 @@ void MagGeoBuilderFromDDD::volumeHandle::buildTruncTubs(const DDExpandedView &fv // will be its center, while this is not true for the largest phi face. buildPhiZSurf(startPhi, deltaPhi, zhalf, rCentr); - if (MagGeoBuilderFromDDD::debug) { - cout << "pos_Rplane " << pos_Rplane << " " << pos_Rplane.perp() << " " << pos_Rplane.phi() << endl - << "rot_R " << surfaces[plane_side]->toGlobal(LocalVector(0., 0., 1.)) << " phi " - << surfaces[plane_side]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl - << "cyl radius " << rCyl << endl; + if (debug) { + LogDebug("magneticfield::volumeHandle") + << "pos_Rplane " << pos_Rplane << " " << pos_Rplane.perp() << " " << pos_Rplane.phi() << endl + << "rot_R " << surfaces[plane_side]->toGlobal(LocalVector(0., 0., 1.)) << " phi " + << surfaces[plane_side]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl + << "cyl radius " << rCyl << endl; // // Check ordering. if ((pos_Rplane.perp() < rCyl) != cutInside) { - cout << "*** WARNING: pos_outer < pos_inner " << endl; + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pos_outer < pos_inner " << endl; } } - // Save volume boundaries theRMin = rIn; theRMax = rOut; diff --git a/MagneticField/GeomBuilder/src/buildTubs.icc b/MagneticField/GeomBuilder/src/buildTubs.icc index b376d1baf3cd9..a06d15a67e3d2 100644 --- a/MagneticField/GeomBuilder/src/buildTubs.icc +++ b/MagneticField/GeomBuilder/src/buildTubs.icc @@ -1,21 +1,18 @@ -void MagGeoBuilderFromDDD::volumeHandle::buildTubs(const DDExpandedView &fv) { - if (MagGeoBuilderFromDDD::debug) - cout << "Building tubs surfaces...: " << endl; +void volumeHandle::buildTubs() { + LogDebug("magneticfield::volumeHandle") << "Building tubs surfaces...: " << endl; DDTubs tubs(solid); - double zhalf = tubs.zhalf() / cm; - double rIn = tubs.rIn() / cm; - double rOut = tubs.rOut() / cm; + // Old DD needs mm to cm conversion, but DD4hep needs no conversion. + // convertUnits should be defined appropriately. + double zhalf = convertUnits(tubs.zhalf()); + double rIn = convertUnits(tubs.rIn()); + double rOut = convertUnits(tubs.rOut()); double startPhi = tubs.startPhi(); double deltaPhi = tubs.deltaPhi(); - if (MagGeoBuilderFromDDD::debug) - cout << "zhalf " << zhalf << endl - << "rIn " << rIn << endl - << "rOut " << rOut << endl - << "startPhi " << startPhi << endl - << "deltaPhi " << deltaPhi << endl; + LogDebug("magneticfield::volumeHandle") << " zhalf " << zhalf << "\n rIn " << rIn << "\n rOut " << rOut + << "\n startPhi " << startPhi << "\n deltaPhi " << deltaPhi; // recalculate center: (for a DDTubs, DDD gives 0,0,Z) double rCentr = (rIn + rOut) / 2.; @@ -36,13 +33,12 @@ void MagGeoBuilderFromDDD::volumeHandle::buildTubs(const DDExpandedView &fv) { buildPhiZSurf(startPhi, deltaPhi, zhalf, rCentr); // Check ordering. - if (MagGeoBuilderFromDDD::debug) { - if (dynamic_cast(&(*surfaces[outer]))->radius() < - dynamic_cast(&(*surfaces[inner]))->radius()) { - cout << "*** WARNING: pos_outer < pos_inner " << endl; + if (debug) { + if (dynamic_cast(&(*surfaces[outer]))->radius() < + dynamic_cast(&(*surfaces[inner]))->radius()) { + LogDebug("magneticfield::volumeHandle") << "*** WARNING: pos_outer < pos_inner " << endl; } } - // Save volume boundaries theRMin = rIn; theRMax = rOut; diff --git a/MagneticField/GeomBuilder/src/eLayer.cc b/MagneticField/GeomBuilder/src/eLayer.cc index e0e5e8d767180..775497275e0dc 100644 --- a/MagneticField/GeomBuilder/src/eLayer.cc +++ b/MagneticField/GeomBuilder/src/eLayer.cc @@ -14,10 +14,10 @@ using namespace SurfaceOrientation; using namespace std; +using namespace magneticfield; //The ctor is in charge of finding sectors inside the layer. -MagGeoBuilderFromDDD::eLayer::eLayer(handles::const_iterator begin, handles::const_iterator end) - : theVolumes(begin, end), mlayer(nullptr) { +eLayer::eLayer(handles::const_iterator begin, handles::const_iterator end) : theVolumes(begin, end), mlayer(nullptr) { // bool debug=MagGeoBuilderFromDDD::debug; // Sort in R @@ -29,8 +29,6 @@ MagGeoBuilderFromDDD::eLayer::eLayer(handles::const_iterator begin, handles::con // } } -MagGeoBuilderFromDDD::eLayer::~eLayer() {} - // double MagGeoBuilderFromDDD::eLayer::minR() const { // // ASSUMPTION: a layer is only 1 volume thick (by construction). // return theVolumes.front()->minR(); @@ -41,7 +39,7 @@ MagGeoBuilderFromDDD::eLayer::~eLayer() {} // return theVolumes.front()->maxR(); // } -MagELayer* MagGeoBuilderFromDDD::eLayer::buildMagELayer() const { +MagELayer* eLayer::buildMagELayer() const { if (mlayer == nullptr) { //FIXME not guaranteed that all volumes in layer have the same zmin // and zmax! diff --git a/MagneticField/GeomBuilder/src/eLayer.h b/MagneticField/GeomBuilder/src/eLayer.h index 726a7ad2cd2d9..1f2196424eb5d 100644 --- a/MagneticField/GeomBuilder/src/eLayer.h +++ b/MagneticField/GeomBuilder/src/eLayer.h @@ -7,28 +7,29 @@ * \author N. Amapane - INFN Torino */ -#include "MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.h" -#include "MagneticField/GeomBuilder/src/volumeHandle.h" #include "MagneticField/GeomBuilder/src/bSector.h" class MagELayer; -class MagGeoBuilderFromDDD::eLayer { -public: - /// Constructor from list of volumes - eLayer(handles::const_iterator begin, handles::const_iterator end); +namespace magneticfield { + class eLayer { + public: + /// Constructor from list of volumes + eLayer(handles::const_iterator begin, handles::const_iterator end); - /// Destructor - ~eLayer(); + /// Destructor + ~eLayer() = default; - // /// Return the list of all volumes. - // const handles & volumes() const {return theVolumes;} + // /// Return the list of all volumes. + // const handles & volumes() const {return theVolumes;} - /// Construct the MagELayer upon request. - MagELayer* buildMagELayer() const; + /// Construct the MagELayer upon request. + MagELayer* buildMagELayer() const; + + private: + handles theVolumes; // pointer to all volumes in this layer + mutable MagELayer* mlayer; + }; +} // namespace magneticfield -private: - handles theVolumes; // pointer to all volumes in this layer - mutable MagELayer* mlayer; -}; #endif diff --git a/MagneticField/GeomBuilder/src/eSector.cc b/MagneticField/GeomBuilder/src/eSector.cc index 1572d27f3851f..60744bf3114b7 100644 --- a/MagneticField/GeomBuilder/src/eSector.cc +++ b/MagneticField/GeomBuilder/src/eSector.cc @@ -7,6 +7,7 @@ */ #include "MagneticField/GeomBuilder/src/eSector.h" +#include "MagneticField/GeomBuilder/src/printUniqueNames.h" #include "Utilities/BinningTools/interface/ClusterizingHistogram.h" #include "MagneticField/Layers/interface/MagESector.h" #include "MagneticField/Layers/interface/MagVerbosity.h" @@ -16,10 +17,11 @@ using namespace SurfaceOrientation; using namespace std; +using namespace magneticfield; // The ctor is in charge of finding layers inside the sector. -MagGeoBuilderFromDDD::eSector::eSector(handles::const_iterator begin, handles::const_iterator end) - : theVolumes(begin, end), msector(nullptr) { +eSector::eSector(handles::const_iterator begin, handles::const_iterator end, bool debugFlag) + : theVolumes(begin, end), msector(nullptr), debug(debugFlag) { //FIXME!!! //precomputed_value_sort(theVolumes.begin(), theVolumes.end(), ExtractAbsZ()); precomputed_value_sort(theVolumes.begin(), theVolumes.end(), ExtractZ()); @@ -30,7 +32,7 @@ MagGeoBuilderFromDDD::eSector::eSector(handles::const_iterator begin, handles::c float zmax = theVolumes.back()->center().z() + resolution; ClusterizingHistogram hisZ(int((zmax - zmin) / resolution) + 1, zmin, zmax); - if (MagGeoBuilderFromDDD::debug) + if (debug) cout << " Z layers: " << zmin << " " << zmax << endl; handles::const_iterator first = theVolumes.begin(); @@ -41,7 +43,7 @@ MagGeoBuilderFromDDD::eSector::eSector(handles::const_iterator begin, handles::c } vector zClust = hisZ.clusterize(resolution); - if (MagGeoBuilderFromDDD::debug) + if (debug) cout << " Found " << zClust.size() << " clusters in Z, " << " layers: " << endl; @@ -52,18 +54,18 @@ MagGeoBuilderFromDDD::eSector::eSector(handles::const_iterator begin, handles::c float zSepar = (zClust[i] + zClust[i + 1]) / 2.f; while ((*separ)->center().z() < zSepar) ++separ; - if (MagGeoBuilderFromDDD::debug) { + if (debug) { cout << " Layer at: " << zClust[i] << " elements: " << separ - layStart << " unique volumes: "; - volumeHandle::printUniqueNames(layStart, separ); + printUniqueNames(layStart, separ); } layers.push_back(eLayer(layStart, separ)); layStart = separ; } { - if (MagGeoBuilderFromDDD::debug) { + if (debug) { cout << " Layer at: " << zClust.back() << " elements: " << last - separ << " unique volumes: "; - volumeHandle::printUniqueNames(separ, last); + printUniqueNames(separ, last); } layers.push_back(eLayer(separ, last)); } @@ -71,9 +73,7 @@ MagGeoBuilderFromDDD::eSector::eSector(handles::const_iterator begin, handles::c // FIXME: Check that all layers have the same dz?. } -MagGeoBuilderFromDDD::eSector::~eSector() {} - -MagESector* MagGeoBuilderFromDDD::eSector::buildMagESector() const { +MagESector* eSector::buildMagESector() const { if (msector == nullptr) { vector mLayers; for (vector::const_iterator lay = layers.begin(); lay != layers.end(); ++lay) { diff --git a/MagneticField/GeomBuilder/src/eSector.h b/MagneticField/GeomBuilder/src/eSector.h index d99fe59388ae2..e89ded20878f8 100644 --- a/MagneticField/GeomBuilder/src/eSector.h +++ b/MagneticField/GeomBuilder/src/eSector.h @@ -8,29 +8,31 @@ * \author N. Amapane - INFN Torino */ -#include "MagneticField/GeomBuilder/src/MagGeoBuilderFromDDD.h" -#include "MagneticField/GeomBuilder/src/volumeHandle.h" #include "MagneticField/GeomBuilder/src/eLayer.h" class MagESector; -class MagGeoBuilderFromDDD::eSector { -public: - /// Constructor from list of volumes - eSector(handles::const_iterator begin, handles::const_iterator end); +namespace magneticfield { + class eSector { + public: + /// Constructor from list of volumes + eSector(handles::const_iterator begin, handles::const_iterator end, bool debugFlag = false); - /// Destructor - ~eSector(); + /// Destructor + ~eSector() = default; - // /// Return all volumes in this sector - // const handles & getVolumes() const {return volumes;} + // /// Return all volumes in this sector + // const handles & getVolumes() const {return volumes;} - /// Construct the MagESector upon request. - MagESector* buildMagESector() const; + /// Construct the MagESector upon request. + MagESector* buildMagESector() const; + + private: + std::vector layers; // the layers in this sectors + handles theVolumes; // pointers to all volumes in the sector + mutable MagESector* msector; + const bool debug; + }; +} // namespace magneticfield -private: - std::vector layers; // the layers in this sectors - handles theVolumes; // pointers to all volumes in the sector - mutable MagESector* msector; -}; #endif diff --git a/MagneticField/GeomBuilder/src/printUniqueNames.cc b/MagneticField/GeomBuilder/src/printUniqueNames.cc new file mode 100644 index 0000000000000..1beb2aecf9048 --- /dev/null +++ b/MagneticField/GeomBuilder/src/printUniqueNames.cc @@ -0,0 +1,29 @@ +#include "MagneticField/GeomBuilder/src/printUniqueNames.h" + +#include +#include +#include + +using namespace std; + +void magneticfield::printUniqueNames(handles::const_iterator begin, handles::const_iterator end, bool uniq) { + std::vector names; + for (handles::const_iterator i = begin; i != end; ++i) { + if (uniq) + names.push_back((*i)->name); + else + names.push_back((*i)->name + ":" + std::to_string((*i)->copyno)); + } + + sort(names.begin(), names.end()); + if (uniq) { + std::vector::iterator i = unique(names.begin(), names.end()); + int nvols = int(i - names.begin()); + cout << nvols << " "; + copy(names.begin(), i, ostream_iterator(cout, " ")); + } else { + cout << names.size() << " "; + copy(names.begin(), names.end(), ostream_iterator(cout, " ")); + } + cout << endl; +} diff --git a/MagneticField/GeomBuilder/src/printUniqueNames.h b/MagneticField/GeomBuilder/src/printUniqueNames.h new file mode 100644 index 0000000000000..0abfe31ce949f --- /dev/null +++ b/MagneticField/GeomBuilder/src/printUniqueNames.h @@ -0,0 +1,12 @@ +#ifndef MagneticField_GeomBuilder_printUniqueNames_h +#define MagneticField_GeomBuilder_printUniqueNames_h + +#include "MagneticField/GeomBuilder/src/BaseVolumeHandle.h" + +namespace magneticfield { + + /// Just for debugging... + void printUniqueNames(handles::const_iterator begin, handles::const_iterator end, bool uniq = true); +} // namespace magneticfield + +#endif diff --git a/MagneticField/GeomBuilder/src/volumeHandle.cc b/MagneticField/GeomBuilder/src/volumeHandle.cc index e1f0ad33e7ed0..24eb0f9ed1fef 100644 --- a/MagneticField/GeomBuilder/src/volumeHandle.cc +++ b/MagneticField/GeomBuilder/src/volumeHandle.cc @@ -31,21 +31,14 @@ using namespace SurfaceOrientation; using namespace std; -MagGeoBuilderFromDDD::volumeHandle::~volumeHandle() { delete refPlane; } - -MagGeoBuilderFromDDD::volumeHandle::volumeHandle(const DDExpandedView &fv, bool expand2Pi) - : name(fv.logicalPart().name().name()), - copyno(fv.copyno()), - magVolume(nullptr), - masterSector(1), - theRN(0.), - theRMin(0.), - theRMax(0.), - refPlane(nullptr), - solid(fv.logicalPart().solid()), - center_(GlobalPoint(fv.translation().x() / cm, fv.translation().y() / cm, fv.translation().z() / cm)), - expand(expand2Pi), - isIronFlag(false) { +MagGeoBuilderFromDDD::volumeHandle::volumeHandle(const DDExpandedView &fv, bool expand2Pi, bool debugVal) + : magneticfield::BaseVolumeHandle(debugVal) { + name = fv.logicalPart().name().name(); + copyno = fv.copyno(); + solid = fv.logicalPart().solid(); + center_ = GlobalPoint(fv.translation().x() / cm, fv.translation().y() / cm, fv.translation().z() / cm); + expand = expand2Pi; + // ASSUMPTION: volume names ends with "_NUM" where NUM is the volume number string volName = name; volName.erase(0, volName.rfind('_') + 1); @@ -55,24 +48,24 @@ MagGeoBuilderFromDDD::volumeHandle::volumeHandle(const DDExpandedView &fv, bool isAssigned[i] = false; } - if (MagGeoBuilderFromDDD::debug) { + if (debug) { cout.precision(7); } referencePlane(fv); if (solid.shape() == DDSolidShape::ddbox) { - buildBox(fv); + buildBox(); } else if (solid.shape() == DDSolidShape::ddtrap) { - buildTrap(fv); + buildTrap(); } else if (solid.shape() == DDSolidShape::ddcons) { - buildCons(fv); + buildCons(); } else if (solid.shape() == DDSolidShape::ddtubs) { - buildTubs(fv); + buildTubs(); } else if (solid.shape() == DDSolidShape::ddpseudotrap) { - buildPseudoTrap(fv); + buildPseudoTrap(); } else if (solid.shape() == DDSolidShape::ddtrunctubs) { - buildTruncTubs(fv); + buildTruncTubs(); } else { cout << "volumeHandle ctor: Unexpected solid: " << DDSolidShapesName::name(solid.shape()) << endl; } @@ -132,7 +125,7 @@ MagGeoBuilderFromDDD::volumeHandle::volumeHandle(const DDExpandedView &fv, bool if (fv.logicalPart().material().name().name() == "Iron") isIronFlag = true; - if (MagGeoBuilderFromDDD::debug) { + if (debug) { cout << " RMin = " << theRMin << endl; cout << " RMax = " << theRMax << endl; @@ -153,8 +146,6 @@ MagGeoBuilderFromDDD::volumeHandle::volumeHandle(const DDExpandedView &fv, bool } } -const Surface::GlobalPoint &MagGeoBuilderFromDDD::volumeHandle::center() const { return center_; } - void MagGeoBuilderFromDDD::volumeHandle::referencePlane(const DDExpandedView &fv) { // The refPlane is the "main plane" for the solid. It corresponds to the // x,y plane in the DDD local frame, and defines a frame where the local @@ -190,7 +181,7 @@ void MagGeoBuilderFromDDD::volumeHandle::referencePlane(const DDExpandedView &fv // The reference plane rotation DD3Vector x, y, z; fv.rotation().GetComponents(x, y, z); - if (MagGeoBuilderFromDDD::debug) { + if (debug) { if (x.Cross(y).Dot(z) < 0.5) { cout << "*** WARNING: Rotation is not RH " << endl; } @@ -210,7 +201,7 @@ void MagGeoBuilderFromDDD::volumeHandle::referencePlane(const DDExpandedView &fv refPlane = new GloballyPositioned(posResult, rotResult); // Check correct orientation - if (MagGeoBuilderFromDDD::debug) { + if (debug) { cout << "Refplane pos " << refPlane->position() << endl; // See comments above for the conventions for orientation. @@ -228,182 +219,6 @@ void MagGeoBuilderFromDDD::volumeHandle::referencePlane(const DDExpandedView &fv } } -void MagGeoBuilderFromDDD::volumeHandle::buildPhiZSurf(double startPhi, double deltaPhi, double zhalf, double rCentr) { - // This is 100% equal for cons and tubs!!! - - GlobalVector planeXAxis = refPlane->toGlobal(LocalVector(1, 0, 0)); - GlobalVector planeYAxis = refPlane->toGlobal(LocalVector(0, 1, 0)); - GlobalVector planeZAxis = refPlane->toGlobal(LocalVector(0, 0, 1)); - - // Local Y axis of the faces at +-phi. - GlobalVector y_phiplus = refPlane->toGlobal(LocalVector(cos(startPhi + deltaPhi), sin(startPhi + deltaPhi), 0.)); - GlobalVector y_phiminus = refPlane->toGlobal(LocalVector(cos(startPhi), sin(startPhi), 0.)); - - Surface::RotationType rot_Z(planeXAxis, planeYAxis); - Surface::RotationType rot_phiplus(planeZAxis, y_phiplus); - Surface::RotationType rot_phiminus(planeZAxis, y_phiminus); - - GlobalPoint pos_zplus(center_.x(), center_.y(), center_.z() + zhalf); - GlobalPoint pos_zminus(center_.x(), center_.y(), center_.z() - zhalf); - // BEWARE: in this case, the origin for phiplus,phiminus surfaces is - // at radius R and not on a plane passing by center_ orthogonal to the radius. - GlobalPoint pos_phiplus( - refPlane->toGlobal(LocalPoint(rCentr * cos(startPhi + deltaPhi), rCentr * sin(startPhi + deltaPhi), 0.))); - GlobalPoint pos_phiminus(refPlane->toGlobal(LocalPoint(rCentr * cos(startPhi), rCentr * sin(startPhi), 0.))); - surfaces[zplus] = new Plane(pos_zplus, rot_Z); - surfaces[zminus] = new Plane(pos_zminus, rot_Z); - surfaces[phiplus] = new Plane(pos_phiplus, rot_phiplus); - surfaces[phiminus] = new Plane(pos_phiminus, rot_phiminus); - - if (MagGeoBuilderFromDDD::debug) { - cout << "Actual Center at: " << center_ << " R " << center_.perp() << " phi " << center_.phi() << endl; - cout << "RN " << theRN << endl; - - cout << "pos_zplus " << pos_zplus << " " << pos_zplus.perp() << " " << pos_zplus.phi() << endl - << "pos_zminus " << pos_zminus << " " << pos_zminus.perp() << " " << pos_zminus.phi() << endl - << "pos_phiplus " << pos_phiplus << " " << pos_phiplus.perp() << " " << pos_phiplus.phi() << endl - << "pos_phiminus " << pos_phiminus << " " << pos_phiminus.perp() << " " << pos_phiminus.phi() << endl; - - cout << "y_phiplus " << y_phiplus << endl; - cout << "y_phiminus " << y_phiminus << endl; - - cout << "rot_Z " << surfaces[zplus]->toGlobal(LocalVector(0., 0., 1.)) << endl - << "rot_phi+ " << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " - << surfaces[phiplus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl - << "rot_phi- " << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)) << " phi " - << surfaces[phiminus]->toGlobal(LocalVector(0., 0., 1.)).phi() << endl; - } - - // // Check ordering. - if (MagGeoBuilderFromDDD::debug) { - if (pos_zplus.z() < pos_zminus.z()) { - cout << "*** WARNING: pos_zplus < pos_zminus " << endl; - } - if (Geom::Phi(pos_phiplus.phi() - pos_phiminus.phi()) < 0.) { - cout << "*** WARNING: pos_phiplus < pos_phiminus " << endl; - } - } -} - -bool MagGeoBuilderFromDDD::volumeHandle::sameSurface(const Surface &s1, Sides which_side, float tolerance) { - //Check for null comparison - if (&s1 == (surfaces[which_side]).get()) { - if (MagGeoBuilderFromDDD::debug) - cout << " sameSurface: OK (same ptr)" << endl; - return true; - } - - const float maxtilt = 0.999; - - const Surface &s2 = *(surfaces[which_side]); - // Try with a plane. - const Plane *p1 = dynamic_cast(&s1); - if (p1 != nullptr) { - const Plane *p2 = dynamic_cast(&s2); - if (p2 == nullptr) { - if (MagGeoBuilderFromDDD::debug) - cout << " sameSurface: different types" << endl; - return false; - } - - if ((fabs(p1->normalVector().dot(p2->normalVector())) > maxtilt) && - (fabs((p1->toLocal(p2->position())).z()) < tolerance)) { - if (MagGeoBuilderFromDDD::debug) - cout << " sameSurface: OK " << fabs(p1->normalVector().dot(p2->normalVector())) << " " - << fabs((p1->toLocal(p2->position())).z()) << endl; - return true; - } else { - if (MagGeoBuilderFromDDD::debug) - cout << " sameSurface: not the same: " << p1->normalVector() << p1->position() << endl - << " " << p2->normalVector() << p2->position() << endl - << fabs(p1->normalVector().dot(p2->normalVector())) << " " << (p1->toLocal(p2->position())).z() << endl; - return false; - } - } - - // Try with a cylinder. - const Cylinder *cy1 = dynamic_cast(&s1); - if (cy1 != nullptr) { - const Cylinder *cy2 = dynamic_cast(&s2); - if (cy2 == nullptr) { - if (MagGeoBuilderFromDDD::debug) - cout << " sameSurface: different types" << endl; - return false; - } - // Assume axis is the same! - if (fabs(cy1->radius() - cy2->radius()) < tolerance) { - return true; - } else { - return false; - } - } - - // Try with a cone. - const Cone *co1 = dynamic_cast(&s1); - if (co1 != nullptr) { - const Cone *co2 = dynamic_cast(&s2); - if (co2 == nullptr) { - if (MagGeoBuilderFromDDD::debug) - cout << " sameSurface: different types" << endl; - return false; - } - // FIXME - if (fabs(co1->openingAngle() - co2->openingAngle()) < maxtilt && - (co1->vertex() - co2->vertex()).mag() < tolerance) { - return true; - } else { - return false; - } - } - - if (MagGeoBuilderFromDDD::debug) - cout << " sameSurface: unknown surfaces..." << endl; - return false; -} - -bool MagGeoBuilderFromDDD::volumeHandle::setSurface(const Surface &s1, Sides which_side) { - //Check for null assignment - if (&s1 == (surfaces[which_side]).get()) { - isAssigned[which_side] = true; - return true; - } - - if (!sameSurface(s1, which_side)) { - cout << "***ERROR: setSurface: trying to assign a surface that does not match destination surface. Skipping." - << endl; - const Surface &s2 = *(surfaces[which_side]); - //FIXME: Just planes for the time being!!! - const Plane *p1 = dynamic_cast(&s1); - const Plane *p2 = dynamic_cast(&s2); - if (p1 != nullptr && p2 != nullptr) - cout << p1->normalVector() << p1->position() << endl << p2->normalVector() << p2->position() << endl; - return false; - } - - if (isAssigned[which_side]) { - if (&s1 != (surfaces[which_side]).get()) { - cout << "*** WARNING volumeHandle::setSurface: trying to reassign a surface to a different surface instance" - << endl; - return false; - } - } else { - surfaces[which_side] = &s1; - isAssigned[which_side] = true; - if (MagGeoBuilderFromDDD::debug) - cout << " Volume " << name << " # " << copyno << " Assigned: " << (int)which_side << endl; - return true; - } - - return false; // let the compiler be happy -} - -const Surface &MagGeoBuilderFromDDD::volumeHandle::surface(Sides which_side) const { return *(surfaces[which_side]); } - -const Surface &MagGeoBuilderFromDDD::volumeHandle::surface(int which_side) const { - assert(which_side >= 0 && which_side < 6); - return *(surfaces[which_side]); -} - std::vector MagGeoBuilderFromDDD::volumeHandle::sides() const { std::vector result; for (int i = 0; i < 6; ++i) { @@ -422,28 +237,15 @@ std::vector MagGeoBuilderFromDDD::volumeHandle::sides() const { return result; } -void MagGeoBuilderFromDDD::volumeHandle::printUniqueNames(handles::const_iterator begin, - handles::const_iterator end, - bool uniq) { - std::vector names; - for (handles::const_iterator i = begin; i != end; ++i) { - if (uniq) - names.push_back((*i)->name); - else - names.push_back((*i)->name + ":" + std::to_string((*i)->copyno)); - } +#include "DataFormats/Math/interface/GeantUnits.h" - sort(names.begin(), names.end()); - if (uniq) { - std::vector::iterator i = unique(names.begin(), names.end()); - int nvols = int(i - names.begin()); - cout << nvols << " "; - copy(names.begin(), i, ostream_iterator(cout, " ")); - } else { - cout << names.size() << " "; - copy(names.begin(), names.end(), ostream_iterator(cout, " ")); - } - cout << endl; +using volumeHandle = MagGeoBuilderFromDDD::volumeHandle; + +// Old DD returns lengths in mm, but CMS code uses cm +template +inline constexpr NumType convertUnits(NumType millimeters) // Millimeters -> centimeters +{ + return (geant_units::operators::convertMmToCm(millimeters)); } #include "MagneticField/GeomBuilder/src/buildBox.icc" diff --git a/MagneticField/GeomBuilder/src/volumeHandle.h b/MagneticField/GeomBuilder/src/volumeHandle.h index 440d84c81b156..2518acd9ecfc9 100644 --- a/MagneticField/GeomBuilder/src/volumeHandle.h +++ b/MagneticField/GeomBuilder/src/volumeHandle.h @@ -17,211 +17,44 @@ #include "DataFormats/GeometrySurface/interface/Surface.h" //#include "ClassReuse/SurfaceGeometry/interface/BoundPlane.h" #include "MagneticField/VolumeGeometry/interface/VolumeSide.h" +#include "MagneticField/GeomBuilder/src/BaseVolumeHandle.h" class DDExpandedView; class MagVolume6Faces; -class MagGeoBuilderFromDDD::volumeHandle { +class MagGeoBuilderFromDDD::volumeHandle : public magneticfield::BaseVolumeHandle { public: - typedef Surface::GlobalPoint GlobalPoint; - typedef Surface::LocalPoint LocalPoint; - typedef Surface::LocalVector LocalVector; - typedef SurfaceOrientation::GlobalFace Sides; - volumeHandle(const DDExpandedView& fv, bool expand2Pi = false); - ~volumeHandle(); + volumeHandle(const DDExpandedView& fv, bool expand2Pi = false, bool debugVal = false); - /// Return the center of the volume - const GlobalPoint& center() const; - /// Distance of (x,y) plane from origin - const double RN() const { return theRN; } - /// Get the current surface on specified side. - const Surface& surface(int which_side) const; - const Surface& surface(Sides which_side) const; - /// Find out if two surfaces are the same physical surface - bool sameSurface(const Surface& s1, Sides which_side, float tolerance = 0.01); - /// Assign a shared surface perorming sanity checks. - bool setSurface(const Surface& s1, Sides which_side); - /// if the specified surface has been matched. - bool isPlaneMatched(int which_side) const { return isAssigned[which_side]; } - - int references(int which_side) const { // FIXME! - /* return surfaces[which_side]->references(); */ - return 0; - } - - /// Name of the volume - std::string name; - /// Name of magnetic field table file - std::string magFile; - /// volume number - unsigned short volumeno; - /// copy number - unsigned short copyno; - - /// Just for debugging... - static void printUniqueNames(handles::const_iterator begin, handles::const_iterator end, bool uniq = true); - - // Phi ranges: Used by: LessDPhiMax; bSector; bSlab::[min|max]Phi(); - // MagBSector, MagBRod - - /// Minimum value of phi covered by the volume - // FIXME: actually returns phi of the point on median plane of the -phi - // surface, except for trapezoids where the absoulte min has been implemented - Geom::Phi minPhi() const { return thePhiMin; } - /// Maximum value of phi covered by the volume - // FIXME: actually returns phi of the point on median plane of the +phi - // surface - Geom::Phi maxPhi() const { return surface(SurfaceOrientation::phiplus).position().phi(); } - - /// Z limits. - // ASSUMPTION: Computed on median Z plane, but this is not a problem since - // all Z planes are orthogonal to the beam line in the current geometry. - double minZ() const { return surface(SurfaceOrientation::zminus).position().z(); } - double maxZ() const { return surface(SurfaceOrientation::zplus).position().z(); } - - /// Minimum R for any point within the volume - double minR() const { return theRMin; } - - /// FIXME: currently returns max RN (to be fixed?). Used by: bLayer::maxR() - // double maxR() const {return theRMax;} - - /// Position and rotation - const GloballyPositioned* placement() const { return refPlane; } - - /// Shape of the solid - DDSolidShape shape() const { return solid.shape(); } - - /// The surfaces and they orientation, as required to build a MagVolume. - std::vector sides() const; - - /// Pointer to the final MagVolume (must be set from outside) - MagVolume6Faces* magVolume; - - bool toExpand() const { return expand; } - - /// Temporary hack to pass information on material. Will eventually be replaced! - bool isIron() const { return isIronFlag; } - - /// The sector for which an interpolator for this class of volumes should be built - int masterSector; - -private: // Disallow Default/copy ctor & assignment op. // (we want to handle only pointers!!!) volumeHandle(const volumeHandle& v) = delete; volumeHandle operator=(const volumeHandle& v) = delete; - // The volume's six surfaces. - RCPS surfaces[6]; - // If the corresponding surface has been assigned to a shared surface. - bool isAssigned[6]; + DDSolidShape shape() const override { return solid.shape(); } + /// The surfaces and they orientation, as required to build a MagVolume. + std::vector sides() const override; + +private: // initialise the refPlane void referencePlane(const DDExpandedView& fv); + // Build the surfaces for a box - void buildBox(const DDExpandedView& fv); + void buildBox(); // Build the surfaces for a trapezoid - void buildTrap(const DDExpandedView& fv); + void buildTrap(); // Build the surfaces for a ddtubs shape - void buildTubs(const DDExpandedView& fv); + void buildTubs(); // Build the surfaces for a ddcons shape - void buildCons(const DDExpandedView& fv); + void buildCons(); // Build the surfaces for a ddpseudotrap shape - void buildPseudoTrap(const DDExpandedView& fv); + void buildPseudoTrap(); // Build the surfaces for a ddtrunctubs shape - void buildTruncTubs(const DDExpandedView& fv); - - // Build phi, z surfaces (common for ddtubs and ddcons) - void buildPhiZSurf(double startPhi, double deltaPhi, double zhalf, double rCentr); - - // Distance from the origin along the normal to the volume's zphi plane. - double theRN; - - // Max and min radius for _any_ point within the volume - // FIXME! - double theRMin; - double theRMax; - Geom::Phi thePhiMin; - - // The refPlane is the "main plane" for the solid. It corresponds to the - // x,y plane in the DDD local frame, and defines a frame where the local - // coordinates are the same as in DDD. - GloballyPositioned* refPlane; + void buildTruncTubs(); // the DDSolid. DDSolid solid; - - // the center of the volume - GlobalPoint center_; - - // Flag this as a master volume out of wich a 2pi volume should be built - // (e.g. central cylinder); this is taken into account by sides(). - bool expand; - - // Temporary hack to keep information on material. Will eventually be replaced! - bool isIronFlag; -}; - -// Extractors for precomputed_value_sort() (safe sorting) - -// To sort volumes in Z -struct MagGeoBuilderFromDDD::ExtractZ { - double operator()(const volumeHandle* v) const { return v->center().z(); } -}; - -// To sort volumes in abs(Z) -struct MagGeoBuilderFromDDD::ExtractAbsZ { - double operator()(const volumeHandle* v) const { return fabs(v->center().z()); } -}; - -// To sort volumes in phi (from -pi to pi). -struct MagGeoBuilderFromDDD::ExtractPhi { - double operator()(const volumeHandle* v) const { - // note that Geom::Phi is implicitly converted to double. - // Periodicity is guaranteed. - return v->center().phi(); - } -}; - -// To sort volumes based on max phi(from -pi to pi). -struct MagGeoBuilderFromDDD::ExtractPhiMax { - double operator()(const volumeHandle* v) const { - // note that Geom::Phi is implicitly converted to double. - // Periodicity is guaranteed. - return v->maxPhi(); - } -}; - -// To sort volumes in R -struct MagGeoBuilderFromDDD::ExtractR { - double operator()(const volumeHandle* v) const { return v->center().perp(); } -}; - -// To sort volumes in RN (distance of (x,y) plane from origin) -struct MagGeoBuilderFromDDD::ExtractRN { - double operator()(const volumeHandle* v) const { return v->RN(); } -}; - -// To sort angles within any range SMALLER THAN PI "counter-clockwise", -// even if the angles cross the pi boundary. -// CAVEAT: // The result is undefined if the input values cover a -// range larger than pi!!! -struct MagGeoBuilderFromDDD::LessDPhi { - bool operator()(double phi1, double phi2) const { - // handle periodicity - return ((Geom::Phi(phi2) - Geom::Phi(phi1)) > 0.); - } -}; - -// Compare the Z of volumes. -// Should be used ONLY for std::max_element and std::min_element -// and NEVER for sorting (use precomputed_value_sort with ExtractZ instead) -struct MagGeoBuilderFromDDD::LessZ { - bool operator()(const volumeHandle* v1, const volumeHandle* v2) const { - if (v1->center().z() < v2->center().z()) - return true; - return false; - } }; #endif diff --git a/MagneticField/GeomBuilder/test/mfwriter.py b/MagneticField/GeomBuilder/test/python/mfwriter.py similarity index 100% rename from MagneticField/GeomBuilder/test/mfwriter.py rename to MagneticField/GeomBuilder/test/python/mfwriter.py diff --git a/MagneticField/GeomBuilder/test/mfxmlwriter.py b/MagneticField/GeomBuilder/test/python/mfxmlwriter.py similarity index 100% rename from MagneticField/GeomBuilder/test/mfxmlwriter.py rename to MagneticField/GeomBuilder/test/python/mfxmlwriter.py diff --git a/MagneticField/GeomBuilder/test/python/testMagGeometry.py b/MagneticField/GeomBuilder/test/python/testMagGeometry.py new file mode 100644 index 0000000000000..199c9c1e5eb2b --- /dev/null +++ b/MagneticField/GeomBuilder/test/python/testMagGeometry.py @@ -0,0 +1,93 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("MagneticFieldTest") + +process.source = cms.Source("EmptySource") +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(2) + ) + +process.load("FWCore.MessageLogger.MessageLogger_cfi") + +process.MessageLogger = cms.Service("MessageLogger", + debugModules = cms.untracked.vstring( '*' ), + destinations = cms.untracked.vstring('cout'), + categories = cms.untracked.vstring ( '*' ), + cout = cms.untracked.PSet( + noLineBreaks = cms.untracked.bool(True), + INFO = cms.untracked.PSet (limit = cms.untracked.int32(-1)), + DEBUG = cms.untracked.PSet (limit = cms.untracked.int32(-1)), + WARNING = cms.untracked.PSet( + limit = cms.untracked.int32(-1) + ), + ERROR = cms.untracked.PSet( + limit = cms.untracked.int32(-1) + ), + threshold = cms.untracked.string('DEBUG'), + default = cms.untracked.PSet (limit = cms.untracked.int32(-1)) + ) +) + + +process.ParametrizedMagneticFieldProducer = cms.ESProducer("ParametrizedMagneticFieldProducer", + version = cms.string('OAE_1103l_071212'), + parameters = cms.PSet( + BValue = cms.string('3_5T') + ), + label = cms.untracked.string('parametrizedField') +) + +process.DDDetectorESProducer = cms.ESSource("DDDetectorESProducer", + confGeomXMLFiles = cms.FileInPath('DetectorDescription/DDCMS/data/cms-mf-geometry.xml'), + appendToDataLabel = cms.string('magfield') + ) + + # version = cms.string('grid_160812_3_8t'), + +process.MagneticFieldESProducer = cms.ESProducer("VolBasedMagFieldESProducerNewDD", + DDDetector = cms.ESInputTag('', 'magfield'), + appendToDataLabel = cms.string(''), + useParametrizedTrackerField = cms.bool(False), + label = cms.untracked.string(''), + attribute = cms.string('magfield'), + value = cms.string('magfield'), + paramLabel = cms.string('parametrizedField'), + version = cms.string('grid_160812_3_5t'), + geometryVersion = cms.int32(160812), + debugBuilder = cms.untracked.bool(True), + cacheLastVolume = cms.untracked.bool(True), + scalingVolumes = cms.vint32(), + scalingFactors = cms.vdouble(), + + gridFiles = cms.VPSet( + cms.PSet( + volumes = cms.string('1001-1010,1012-1027,1030-1033,1036-1041,1044-1049,1052-1057,1060-1063,1066-1071,1074-1077,1080-1097,1102-1129,1138-1402,1415-1416,' + + '2001-2010,2012-2027,2030-2033,2036-2041,2044-2049,2052-2057,2060-2063,2066-2071,2074-2077,2080-2097,2102-2129,2138-2402,2415-2416'), + sectors = cms.string('0') , + master = cms.int32(0), + path = cms.string('s[s]/grid.[v].bin'), + ), + # Replicate sector 1 for volumes outside any detector + cms.PSet( + volumes = cms.string('1011,1028-1029,1034-1035,1042-1043,1050-1051,1058-1059,1064-1065,1072-1073,1078-1079,'+ # volumes extending from R~7.6 m to to R=9 m, + '1098-1101,1130-1137,' + # Forward volumes, ouside CASTOR/HF + '1403-1414,1417-1464,' # Volumes beyond |Z|>17.74 + '2011,2028-2029,2034-2035,2042-2043,2050-2051,2058-2059,2064-2065,2072-2073,2078-2079,'+ + '2098-2101,2130-2137,'+ + '2403-2414,2417-2464'), + sectors = cms.string('0'), + master = cms.int32(1), + path = cms.string('s01/grid.[v].bin'), + ), + ) + ) + +process.DDCompactViewMFESProducer = cms.ESProducer("DDCompactViewMFESProducer", + appendToDataLabel = cms.string('magfield') + ) + +process.test = cms.EDAnalyzer("testMagGeometryAnalyzer", + DDDetector = cms.ESInputTag('', 'magfield') + ) + +process.p = cms.Path(process.test) diff --git a/MagneticField/GeomBuilder/test/python/testMagGeometryOldDD.py b/MagneticField/GeomBuilder/test/python/testMagGeometryOldDD.py new file mode 100644 index 0000000000000..6e04b2cfc74a8 --- /dev/null +++ b/MagneticField/GeomBuilder/test/python/testMagGeometryOldDD.py @@ -0,0 +1,102 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("MagneticFieldTest") + +process.source = cms.Source("EmptySource") +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(2) + ) + +process.load("FWCore.MessageLogger.MessageLogger_cfi") + +process.MessageLogger = cms.Service("MessageLogger", + debugModules = cms.untracked.vstring( '*' ), + destinations = cms.untracked.vstring('cout'), + categories = cms.untracked.vstring ( '*' ), + cout = cms.untracked.PSet( + noLineBreaks = cms.untracked.bool(True), + INFO = cms.untracked.PSet (limit = cms.untracked.int32(-1)), + DEBUG = cms.untracked.PSet (limit = cms.untracked.int32(-1)), + WARNING = cms.untracked.PSet( + limit = cms.untracked.int32(-1) + ), + ERROR = cms.untracked.PSet( + limit = cms.untracked.int32(-1) + ), + threshold = cms.untracked.string('DEBUG'), + default = cms.untracked.PSet (limit = cms.untracked.int32(-1)) + ) +) + +process.magfield = cms.ESSource("XMLIdealGeometryESSource", + geomXMLFiles = cms.vstring('Geometry/CMSCommonData/data/normal/cmsextent.xml', + 'Geometry/CMSCommonData/data/cms.xml', + 'Geometry/CMSCommonData/data/cmsMagneticField.xml', + 'MagneticField/GeomBuilder/data/MagneticFieldVolumes_160812_1.xml', + 'MagneticField/GeomBuilder/data/MagneticFieldVolumes_160812_2.xml', + 'Geometry/CMSCommonData/data/materials.xml'), + rootNodeName = cms.string('cmsMagneticField:MAGF') +) + +# avoid interference with EmptyESSource in uniformMagneticField.cfi +process.es_prefer_magfield = cms.ESPrefer("XMLIdealGeometryESSource","magfield") + + +process.ParametrizedMagneticFieldProducer = cms.ESProducer("ParametrizedMagneticFieldProducer", + version = cms.string('OAE_1103l_071212'), + parameters = cms.PSet( + BValue = cms.string('3_5T') + ), + label = cms.untracked.string('parametrizedField') +) + +process.DDDetectorESProducer = cms.ESSource("DDDetectorESProducer", + confGeomXMLFiles = cms.FileInPath('DetectorDescription/DDCMS/data/cms-mf-geometry.xml'), + appendToDataLabel = cms.string('magfield') + ) + + # version = cms.string('grid_160812_3_8t'), + +process.MagneticFieldESProducer = cms.ESProducer("VolumeBasedMagneticFieldESProducer", + DDDetector = cms.ESInputTag('', 'magfield'), + appendToDataLabel = cms.string(''), + useParametrizedTrackerField = cms.bool(True), + label = cms.untracked.string(''), + attribute = cms.string('magfield'), + value = cms.string('magfield'), + paramLabel = cms.string('parametrizedField'), + version = cms.string('grid_160812_3_5t'), + geometryVersion = cms.int32(160812), + debugBuilder = cms.untracked.bool(True), + cacheLastVolume = cms.untracked.bool(True), + scalingVolumes = cms.vint32(), + scalingFactors = cms.vdouble(), + + gridFiles = cms.VPSet( + cms.PSet( + volumes = cms.string('1001-1010,1012-1027,1030-1033,1036-1041,1044-1049,1052-1057,1060-1063,1066-1071,1074-1077,1080-1097,1102-1129,1138-1402,1415-1416,' + + '2001-2010,2012-2027,2030-2033,2036-2041,2044-2049,2052-2057,2060-2063,2066-2071,2074-2077,2080-2097,2102-2129,2138-2402,2415-2416'), + sectors = cms.string('0') , + master = cms.int32(0), + path = cms.string('s[s]/grid.[v].bin'), + ), + # Replicate sector 1 for volumes outside any detector + cms.PSet( + volumes = cms.string('1011,1028-1029,1034-1035,1042-1043,1050-1051,1058-1059,1064-1065,1072-1073,1078-1079,'+ # volumes extending from R~7.6 m to to R=9 m, + '1098-1101,1130-1137,' + # Forward volumes, ouside CASTOR/HF + '1403-1414,1417-1464,' # Volumes beyond |Z|>17.74 + '2011,2028-2029,2034-2035,2042-2043,2050-2051,2058-2059,2064-2065,2072-2073,2078-2079,'+ + '2098-2101,2130-2137,'+ + '2403-2414,2417-2464'), + sectors = cms.string('0'), + master = cms.int32(1), + path = cms.string('s01/grid.[v].bin'), + ), + ) + ) + +process.test = cms.EDAnalyzer("testMagGeometryAnalyzer", + DDDetector = cms.ESInputTag('', 'magfield') + ) + +process.p = cms.Path(process.test) diff --git a/MagneticField/GeomBuilder/test/stubs/MagGeometryAnalyzer.cc b/MagneticField/GeomBuilder/test/stubs/MagGeometryAnalyzer.cc index 81a94876a1fd6..3d47ae7e55709 100644 --- a/MagneticField/GeomBuilder/test/stubs/MagGeometryAnalyzer.cc +++ b/MagneticField/GeomBuilder/test/stubs/MagGeometryAnalyzer.cc @@ -50,7 +50,8 @@ void testMagGeometryAnalyzer::analyze(const edm::Event& event, const edm::EventS MagGeometryExerciser exe(field); //FIXME: the region to be tested is specified inside. - exe.testFindVolume(10000000); + // exe.testFindVolume(10000000); + exe.testFindVolume(1000); // Test that random points are inside one and only one volume // exe.testInside(100000,0.03); diff --git a/MagneticField/GeomBuilder/test/stubs/MagGeometryExerciser.cc b/MagneticField/GeomBuilder/test/stubs/MagGeometryExerciser.cc index 73023a32bc934..b4ce82424162a 100644 --- a/MagneticField/GeomBuilder/test/stubs/MagGeometryExerciser.cc +++ b/MagneticField/GeomBuilder/test/stubs/MagGeometryExerciser.cc @@ -64,7 +64,7 @@ bool MagGeometryExerciser::testFindVolume(const GlobalPoint& gp) { bool ok = (vol != 0); if (vol == 0) { - cout << "ERROR no volume found! " << gp << " " << gp.z() << " " << gp.perp() + cout << "Test ERROR: No volume found! Global point: " << gp << " , GP Z = " << gp.z() << ", GP perp = " << gp.perp() << " isBarrel: " << theGeometry->inBarrel(gp) << endl; // Try with a linear search diff --git a/MagneticField/Records/interface/IdealMagneticFieldRecord.h b/MagneticField/Records/interface/IdealMagneticFieldRecord.h index d08afe03a1644..1982eee5c6d4e 100644 --- a/MagneticField/Records/interface/IdealMagneticFieldRecord.h +++ b/MagneticField/Records/interface/IdealMagneticFieldRecord.h @@ -6,10 +6,12 @@ #include "CondFormats/DataRecord/interface/RunSummaryRcd.h" #include "CondFormats/DataRecord/interface/MagFieldConfigRcd.h" #include "CondFormats/DataRecord/interface/MFGeometryFileRcd.h" +#include "Geometry/Records/interface/GeometryFileRcd.h" #include "boost/mpl/vector.hpp" -class IdealMagneticFieldRecord : public edm::eventsetup::DependentRecordImplementation< - IdealMagneticFieldRecord, - boost::mpl::vector > {}; +class IdealMagneticFieldRecord + : public edm::eventsetup::DependentRecordImplementation< + IdealMagneticFieldRecord, + boost::mpl::vector > {}; #endif diff --git a/MagneticField/VolumeBasedEngine/interface/MagGeometry.h b/MagneticField/VolumeBasedEngine/interface/MagGeometry.h index 12a22b7e3cab0..1671828522a73 100644 --- a/MagneticField/VolumeBasedEngine/interface/MagGeometry.h +++ b/MagneticField/VolumeBasedEngine/interface/MagGeometry.h @@ -9,7 +9,6 @@ #include "DataFormats/GeometrySurface/interface/BoundPlane.h" #include "MagneticField/Layers/src/MagBinFinders.h" -#include "DetectorDescription/Core/interface/DDCompactView.h" #include #include diff --git a/MagneticField/VolumeBasedEngine/src/MagGeometry.cc b/MagneticField/VolumeBasedEngine/src/MagGeometry.cc index dae67c798a52e..3f11aa0c9a45b 100644 --- a/MagneticField/VolumeBasedEngine/src/MagGeometry.cc +++ b/MagneticField/VolumeBasedEngine/src/MagGeometry.cc @@ -13,10 +13,13 @@ #include "Utilities/BinningTools/interface/PeriodicBinFinderInPhi.h" #include "FWCore/Utilities/interface/isFinite.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" #include "MagneticField/Layers/interface/MagVerbosity.h" #include "FWCore/MessageLogger/interface/MessageLogger.h" +#include + using namespace std; using namespace edm; @@ -63,12 +66,15 @@ MagGeometry::MagGeometry(int geomVersion, //FIXME assume sectors are already sorted in phi //FIXME: PeriodicBinFinderInPhi gets *center* of first bin int nEBins = theESectors.size(); - theEndcapBinFinder = new PeriodicBinFinderInPhi(theESectors.front()->minPhi() + Geom::pi() / nEBins, nEBins); + if (nEBins > 0) + theEndcapBinFinder = new PeriodicBinFinderInPhi(theESectors.front()->minPhi() + Geom::pi() / nEBins, nEBins); } MagGeometry::~MagGeometry() { - delete theBarrelBinFinder; - delete theEndcapBinFinder; + if (theBarrelBinFinder != nullptr) + delete theBarrelBinFinder; + if (theEndcapBinFinder != nullptr) + delete theEndcapBinFinder; for (vector::const_iterator ilay = theBLayers.begin(); ilay != theBLayers.end(); ++ilay) { delete (*ilay); @@ -103,11 +109,16 @@ GlobalVector MagGeometry::fieldInTesla(const GlobalPoint& gp) const { MagVolume const* MagGeometry::findVolume1(const GlobalPoint& gp, double tolerance) const { MagVolume6Faces const* found = nullptr; + int errCnt = 0; if (inBarrel(gp)) { // Barrel for (vector::const_iterator v = theBVolumes.begin(); v != theBVolumes.end(); ++v) { if ((*v) == nullptr) { //FIXME: remove this check - cout << endl << "***ERROR: MagGeometry::findVolume: MagVolume not set" << endl; - continue; + cout << endl << "***ERROR: MagGeometry::findVolume: MagVolume for barrel not set" << endl; + ++errCnt; + if (errCnt < 3) + continue; + else + break; } if ((*v)->inside(gp, tolerance)) { found = (*v); @@ -118,8 +129,12 @@ MagVolume const* MagGeometry::findVolume1(const GlobalPoint& gp, double toleranc } else { // Endcaps for (vector::const_iterator v = theEVolumes.begin(); v != theEVolumes.end(); ++v) { if ((*v) == nullptr) { //FIXME: remove this check - cout << endl << "***ERROR: MagGeometry::findVolume: MagVolume not set" << endl; - continue; + cout << endl << "***ERROR: MagGeometry::findVolume: MagVolume for endcap not set" << endl; + ++errCnt; + if (errCnt < 3) + continue; + else + break; } if ((*v)->inside(gp, tolerance)) { found = (*v); @@ -157,12 +172,15 @@ MagVolume const* MagGeometry::findVolume(const GlobalPoint& gp, double tolerance } else { // Endcaps Geom::Phi phi = gp.phi(); - int bin = theEndcapBinFinder->binIndex(phi); - if (verbose::debugOut) - cout << "Trying endcap sector at phi " << theESectors[bin]->minPhi() << " " << phi << endl; - result = theESectors[bin]->findVolume(gp, tolerance); - if (verbose::debugOut) - cout << "***In guessed esector " << (result == nullptr ? " failed " : " OK ") << endl; + if (theEndcapBinFinder != nullptr && !theESectors.empty()) { + int bin = theEndcapBinFinder->binIndex(phi); + if (verbose::debugOut) + cout << "Trying endcap sector at phi " << theESectors[bin]->minPhi() << " " << phi << endl; + result = theESectors[bin]->findVolume(gp, tolerance); + if (verbose::debugOut) + cout << "***In guessed esector " << (result == nullptr ? " failed " : " OK ") << endl; + } else + edm::LogError("MagGeometry") << "Endcap empty"; } if (result == nullptr && tolerance < 0.0001) {