diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/DistortedMaterialManager.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/DistortedMaterialManager.cxx index d07b148d1f1bc12c9ceabf498ba77b658b0d2afd..c7db36eaa8420a5408ba4b9992754496038a9ee5 100755 --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/DistortedMaterialManager.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/DistortedMaterialManager.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/DistortedMaterialManager.h" @@ -13,39 +13,34 @@ #include "GaudiKernel/Bootstrap.h" namespace InDetDD { + DistortedMaterialManager::DistortedMaterialManager() { + ISvcLocator* svcLocator = Gaudi::svcLocator(); // from Bootstrap -DistortedMaterialManager::DistortedMaterialManager() -{ - - ISvcLocator* svcLocator = Gaudi::svcLocator(); // from Bootstrap + Athena::MsgStreamMember log("ExtraMaterialManager"); + log << MSG::DEBUG << "Initialized InDet Distorted Material Manager" << endmsg; - Athena::MsgStreamMember log("ExtraMaterialManager"); - log << MSG::DEBUG << "Initialized InDet Distorted Material Manager" << endmsg; + StoreGateSvc* detStore; + StatusCode sc; + sc = svcLocator->service("DetectorStore", detStore); + if (sc.isFailure()) log << MSG::FATAL << "Could not locate DetectorStore" << endmsg; - StoreGateSvc * detStore; - StatusCode sc; - sc = svcLocator->service("DetectorStore", detStore ); - if (sc.isFailure()) log << MSG::FATAL << "Could not locate DetectorStore" << endmsg; - - IRDBAccessSvc *rdbSvc; - sc = svcLocator->service("RDBAccessSvc",rdbSvc); - if (sc.isFailure()) log << MSG::FATAL << "Could not locate RDBAccessSvc" << endmsg; + IRDBAccessSvc* rdbSvc; + sc = svcLocator->service("RDBAccessSvc", rdbSvc); + if (sc.isFailure()) log << MSG::FATAL << "Could not locate RDBAccessSvc" << endmsg; - // Get version tag and node for InDet. - DecodeVersionKey versionKey("InnerDetector"); - std::string detectorKey = versionKey.tag(); - std::string detectorNode = versionKey.node(); + // Get version tag and node for InDet. + DecodeVersionKey versionKey("InnerDetector"); + std::string detectorKey = versionKey.tag(); + std::string detectorNode = versionKey.node(); - log << MSG::DEBUG << "Retrieving Record Sets from database ..." << endmsg; - log << MSG::DEBUG << "Key = " << detectorKey << " Node = " << detectorNode << endmsg; + log << MSG::DEBUG << "Retrieving Record Sets from database ..." << endmsg; + log << MSG::DEBUG << "Key = " << detectorKey << " Node = " << detectorNode << endmsg; - m_xMatTable = rdbSvc->getRecordsetPtr("InDetExtraMaterial", detectorKey, detectorNode); - - const StoredMaterialManager * theGeoMaterialManager = 0; - sc = detStore->retrieve(theGeoMaterialManager, "MATERIALS"); - if (sc.isFailure()) log << MSG::FATAL << "Could not locate GeoModel Material manager" << endmsg; - m_materialManager = theGeoMaterialManager; - -} + m_xMatTable = rdbSvc->getRecordsetPtr("InDetExtraMaterial", detectorKey, detectorNode); + const StoredMaterialManager* theGeoMaterialManager = 0; + sc = detStore->retrieve(theGeoMaterialManager, "MATERIALS"); + if (sc.isFailure()) log << MSG::FATAL << "Could not locate GeoModel Material manager" << endmsg; + m_materialManager = theGeoMaterialManager; + } } // end namespace diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ExtraMaterial.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ExtraMaterial.cxx index c5a6e9d9a6784362706455792641123063a76181..87de0f163a37107375ac09a7a35071aeb57a11b6 100755 --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ExtraMaterial.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ExtraMaterial.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/ExtraMaterial.h" #include "InDetGeoModelUtils/GenericTubeMaker.h" @@ -27,62 +27,51 @@ #include <iostream> namespace InDetDD { + ExtraMaterial::ExtraMaterial(IRDBRecordset_ptr xMatTable, const AbsMaterialManager* matManager) + : m_xMatTable(xMatTable), + m_matManager(matManager) + {} + + ExtraMaterial::ExtraMaterial(const DistortedMaterialManager* manager) + : m_xMatTable(manager->extraMaterialTable()), + m_matManager(manager->materialManager()) + {} + + void + ExtraMaterial::add(GeoPhysVol* parent, const std::string& region, double zParent) { + add(parent, 0, region, zParent); + } + void + ExtraMaterial::add(GeoFullPhysVol* parent, const std::string& region, double zParent) { + add(0, parent, region, zParent); + } + void + ExtraMaterial::add(GeoPhysVol* parent, GeoFullPhysVol* fullparent, const std::string& region, double zParent) { + //std::cout << "Adding Extra material for region: " << region << ", zParent = " << zParent << std::endl; -ExtraMaterial::ExtraMaterial(IRDBRecordset_ptr xMatTable, const AbsMaterialManager * matManager) - : m_xMatTable(xMatTable), - m_matManager(matManager) -{} - -ExtraMaterial::ExtraMaterial(const DistortedMaterialManager * manager) - : m_xMatTable(manager->extraMaterialTable()), - m_matManager(manager->materialManager()) -{} - -void -ExtraMaterial::add(GeoPhysVol * parent, const std::string & region, double zParent) -{ - add(parent, 0, region, zParent); -} - -void -ExtraMaterial::add(GeoFullPhysVol * parent, const std::string & region, double zParent) -{ - add(0, parent, region, zParent); -} - -void -ExtraMaterial::add(GeoPhysVol * parent, GeoFullPhysVol * fullparent, const std::string & region, double zParent) -{ + for (unsigned int i = 0; i < m_xMatTable->size(); i++) { + std::ostringstream volnamestr; + volnamestr << "ExtraMaterial" << i; - //std::cout << "Adding Extra material for region: " << region << ", zParent = " << zParent << std::endl; + //std::cout << "In Extra material " << i << std::endl; - for (unsigned int i = 0; i < m_xMatTable->size(); i++) { - std::ostringstream volnamestr; - volnamestr << "ExtraMaterial" << i; - - //std::cout << "In Extra material " << i << std::endl; + if ((*m_xMatTable)[i]->getString("REGION") == region) { + //std::cout << "Extra material Match " << i << std::endl; - if ((*m_xMatTable)[i]->getString("REGION") == region) { + GenericTubeMaker tubeHelper((*m_xMatTable)[i]); + const GeoMaterial* material = m_matManager->getMaterial(tubeHelper.volData().material()); + const GeoShape* shape = tubeHelper.buildShape(); + GeoLogVol* logVol = new GeoLogVol(volnamestr.str(), shape, material); + GeoPhysVol* physVol = new GeoPhysVol(logVol); - //std::cout << "Extra material Match " << i << std::endl; - - GenericTubeMaker tubeHelper((*m_xMatTable)[i]); - const GeoMaterial * material = m_matManager->getMaterial(tubeHelper.volData().material()); - const GeoShape * shape = tubeHelper.buildShape(); - GeoLogVol * logVol = new GeoLogVol(volnamestr.str(), shape, material); - GeoPhysVol * physVol = new GeoPhysVol(logVol); - - if (parent) { - tubeHelper.placeVolume(parent,physVol,zParent); - } else { - tubeHelper.placeVolume(fullparent,physVol,zParent); + if (parent) { + tubeHelper.placeVolume(parent, physVol, zParent); + } else { + tubeHelper.placeVolume(fullparent, physVol, zParent); + } } } } -} - - } // end namespace - diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/GenericTubeMaker.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/GenericTubeMaker.cxx index cbbe27ca1361ddb7eef1b16e5b6d140cb367fa41..920a3f051611a1748435fc7b09f24ccb3e37dba2 100755 --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/GenericTubeMaker.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/GenericTubeMaker.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/GenericTubeMaker.h" #include "InDetGeoModelUtils/TubeVolData.h" @@ -16,154 +16,147 @@ #include "RDBAccessSvc/IRDBRecord.h" namespace InDetDD { - -GenericTubeMaker::GenericTubeMaker(const IRDBRecord * record) - : m_record(record), + GenericTubeMaker::GenericTubeMaker(const IRDBRecord* record) + : m_record(record), m_volData(record) -{} - -std::string -GenericTubeMaker::materialName() const -{ - return m_record->getString("MATERIAL"); -} - -std::string -GenericTubeMaker::name() const -{ - return m_record->getString("NAME"); -} - -const GeoShape * -GenericTubeMaker::buildShape() -{ - const GeoShape * shape = 0; - - switch(m_volData.shape()) { - case TubeVolData::TUBE : - shape = new GeoTube(m_volData.rmin(), m_volData.rmax(), 0.5*m_volData.length()); - break; - case TubeVolData::TUBS : - shape = new GeoTubs(m_volData.rmin(), m_volData.rmax(), 0.5*m_volData.length(), m_volData.phiStart(), m_volData.phiDelta()); - break; - case TubeVolData::CONS : - shape = new GeoCons(m_volData.rmin(), m_volData.rmin2(), m_volData.rmax(), m_volData.rmax2(), - 0.5*m_volData.length(), m_volData.phiStart(), m_volData.phiDelta()); - break; - case TubeVolData::RADIAL : - // This simulates the radial decrease in density. - double zstart = -0.5*m_volData.length(); - GeoPcon * shapeTmp = new GeoPcon(m_volData.phiStart(), m_volData.phiDelta()); - shapeTmp->addPlane(zstart, m_volData.rmin(), m_volData.rmax()); - double radialDelta = (m_volData.rmax() - m_volData.rmin())/m_volData.radialDivisions(); - for (int i = 0; i < m_volData.radialDivisions(); i++){ - double rIntermediate = m_volData.rmax() - i*radialDelta; - double reductionFactor = m_volData.rmin()/rIntermediate; - shapeTmp->addPlane(zstart + reductionFactor*m_volData.length(), m_volData.rmin(), rIntermediate); - } - shapeTmp->addPlane(zstart + m_volData.length(), m_volData.rmin(), m_volData.rmin()); - shape = shapeTmp; - break; + {} + + std::string + GenericTubeMaker::materialName() const { + return m_record->getString("MATERIAL"); } - return shape; -} - -void -GenericTubeMaker::placeVolume(GeoPhysVol * parent, GeoVPhysVol * child, double zParent) -{ - placeVolume(parent, 0, child, zParent); -} - -void -GenericTubeMaker::placeVolume(GeoFullPhysVol * fullparent, GeoVPhysVol * child, double zParent) -{ - placeVolume(0, fullparent, child, zParent); -} - -void -GenericTubeMaker::placeVolTwoSide(GeoPhysVol * parentPos, GeoPhysVol * parentNeg, GeoVPhysVol * child, double zParent) -{ - placeVolTwoSide(parentPos, parentNeg, 0, 0, child, zParent); -} - -void -GenericTubeMaker::placeVolTwoSide(GeoFullPhysVol * fullparentPos, GeoFullPhysVol * fullparentNeg, GeoVPhysVol * child, double zParent) -{ - placeVolTwoSide(0, 0, fullparentPos, fullparentNeg, child, zParent); -} - -void -GenericTubeMaker::placeVolume(GeoPhysVol * parent, GeoFullPhysVol * fullparent, GeoVPhysVol * child, double zParent) -{ - for (int iRepeat = 0; iRepeat < m_volData.nRepeat(); iRepeat++) { - - double phi = m_volData.phiStep() * iRepeat; - - GeoTransform * xform = 0; - double zOffset = m_volData.zMid()-zParent; - if (zOffset != 0 || iRepeat > 0) { - xform = new GeoTransform(HepGeom::TranslateZ3D(zOffset)*HepGeom::RotateZ3D(phi)); - } - - if (parent) { - if (xform) parent->add(xform); - parent->add(child); - } else { - if (xform) fullparent->add(xform); - fullparent->add(child); - } + std::string + GenericTubeMaker::name() const { + return m_record->getString("NAME"); + } - // Place in negative z as well. - if (m_volData.bothZ()) { - GeoTransform * xformNeg = new GeoTransform(HepGeom::RotateY3D(180*CLHEP::deg)*HepGeom::TranslateZ3D(zOffset)*HepGeom::RotateZ3D(phi)); - if (parent) { - parent->add(xformNeg); - parent->add(child); - } else { - fullparent->add(xformNeg); - fullparent->add(child); + const GeoShape* + GenericTubeMaker::buildShape() { + const GeoShape* shape = 0; + + switch (m_volData.shape()) { + case TubeVolData::TUBE: + shape = new GeoTube(m_volData.rmin(), m_volData.rmax(), 0.5 * m_volData.length()); + break; + + case TubeVolData::TUBS: + shape = new GeoTubs(m_volData.rmin(), m_volData.rmax(), 0.5 * m_volData.length(), + m_volData.phiStart(), m_volData.phiDelta()); + break; + + case TubeVolData::CONS: + shape = new GeoCons(m_volData.rmin(), m_volData.rmin2(), m_volData.rmax(), m_volData.rmax2(), + 0.5 * m_volData.length(), m_volData.phiStart(), m_volData.phiDelta()); + break; + + case TubeVolData::RADIAL: + // This simulates the radial decrease in density. + double zstart = -0.5 * m_volData.length(); + GeoPcon* shapeTmp = new GeoPcon(m_volData.phiStart(), m_volData.phiDelta()); + shapeTmp->addPlane(zstart, m_volData.rmin(), m_volData.rmax()); + double radialDelta = (m_volData.rmax() - m_volData.rmin()) / m_volData.radialDivisions(); + for (int i = 0; i < m_volData.radialDivisions(); i++) { + double rIntermediate = m_volData.rmax() - i * radialDelta; + double reductionFactor = m_volData.rmin() / rIntermediate; + shapeTmp->addPlane(zstart + reductionFactor * m_volData.length(), m_volData.rmin(), rIntermediate); } + shapeTmp->addPlane(zstart + m_volData.length(), m_volData.rmin(), m_volData.rmin()); + shape = shapeTmp; + break; } - } // iRepeat loop -} - -void -GenericTubeMaker::placeVolTwoSide(GeoPhysVol * parentPos, GeoPhysVol * parentNeg, - GeoFullPhysVol * fullparentPos, GeoFullPhysVol * fullparentNeg, - GeoVPhysVol * child, double zParent) -{ - for (int iRepeat = 0; iRepeat < m_volData.nRepeat(); iRepeat++) { - double phi = m_volData.phiStep() * iRepeat; - double zOffset = m_volData.zMid()-zParent; - const bool newXform((zOffset != 0) or (iRepeat > 0)); - - if (parentPos) { - if (newXform){ - parentPos->add(new GeoTransform(HepGeom::TranslateZ3D(zOffset)*HepGeom::RotateZ3D(phi))); - } - parentPos->add(child); - } else if(fullparentPos){ - if (newXform){ - fullparentPos->add(new GeoTransform(HepGeom::TranslateZ3D(zOffset)*HepGeom::RotateZ3D(phi))); + return shape; + } + + void + GenericTubeMaker::placeVolume(GeoPhysVol* parent, GeoVPhysVol* child, double zParent) { + placeVolume(parent, 0, child, zParent); + } + + void + GenericTubeMaker::placeVolume(GeoFullPhysVol* fullparent, GeoVPhysVol* child, double zParent) { + placeVolume(0, fullparent, child, zParent); + } + + void + GenericTubeMaker::placeVolTwoSide(GeoPhysVol* parentPos, GeoPhysVol* parentNeg, GeoVPhysVol* child, double zParent) { + placeVolTwoSide(parentPos, parentNeg, 0, 0, child, zParent); + } + + void + GenericTubeMaker::placeVolTwoSide(GeoFullPhysVol* fullparentPos, GeoFullPhysVol* fullparentNeg, GeoVPhysVol* child, + double zParent) { + placeVolTwoSide(0, 0, fullparentPos, fullparentNeg, child, zParent); + } + + void + GenericTubeMaker::placeVolume(GeoPhysVol* parent, GeoFullPhysVol* fullparent, GeoVPhysVol* child, double zParent) { + for (int iRepeat = 0; iRepeat < m_volData.nRepeat(); iRepeat++) { + double phi = m_volData.phiStep() * iRepeat; + + GeoTransform* xform = 0; + double zOffset = m_volData.zMid() - zParent; + if (zOffset != 0 || iRepeat > 0) { + xform = new GeoTransform(HepGeom::TranslateZ3D(zOffset) * HepGeom::RotateZ3D(phi)); } - fullparentPos->add(child); - } - // Place in negative z as well. - if (m_volData.bothZ()) { - GeoTransform * xformNeg = new GeoTransform(HepGeom::RotateY3D(180*CLHEP::deg)*HepGeom::TranslateZ3D(zOffset)*HepGeom::RotateZ3D(phi)); - if (parentNeg) { - parentNeg->add(xformNeg); - parentNeg->add(child); + if (parent) { + if (xform) parent->add(xform); + parent->add(child); } else { - fullparentNeg->add(xformNeg); - fullparentNeg->add(child); + if (xform) fullparent->add(xform); + fullparent->add(child); } - } - } // iRepeat loop -} + // Place in negative z as well. + if (m_volData.bothZ()) { + GeoTransform* xformNeg = new GeoTransform(HepGeom::RotateY3D(180 * CLHEP::deg) * HepGeom::TranslateZ3D( + zOffset) * HepGeom::RotateZ3D(phi)); + if (parent) { + parent->add(xformNeg); + parent->add(child); + } else { + fullparent->add(xformNeg); + fullparent->add(child); + } + } + } // iRepeat loop + } -}// end namespace + void + GenericTubeMaker::placeVolTwoSide(GeoPhysVol* parentPos, GeoPhysVol* parentNeg, + GeoFullPhysVol* fullparentPos, GeoFullPhysVol* fullparentNeg, + GeoVPhysVol* child, double zParent) { + for (int iRepeat = 0; iRepeat < m_volData.nRepeat(); iRepeat++) { + double phi = m_volData.phiStep() * iRepeat; + double zOffset = m_volData.zMid() - zParent; + const bool newXform((zOffset != 0)or(iRepeat > 0)); + + if (parentPos) { + if (newXform) { + parentPos->add(new GeoTransform(HepGeom::TranslateZ3D(zOffset) * HepGeom::RotateZ3D(phi))); + } + parentPos->add(child); + } else if (fullparentPos) { + if (newXform) { + fullparentPos->add(new GeoTransform(HepGeom::TranslateZ3D(zOffset) * HepGeom::RotateZ3D(phi))); + } + fullparentPos->add(child); + } + + // Place in negative z as well. + if (m_volData.bothZ()) { + GeoTransform* xformNeg = new GeoTransform(HepGeom::RotateY3D(180 * CLHEP::deg) * HepGeom::TranslateZ3D( + zOffset) * HepGeom::RotateZ3D(phi)); + if (parentNeg) { + parentNeg->add(xformNeg); + parentNeg->add(child); + } else { + fullparentNeg->add(xformNeg); + fullparentNeg->add(child); + } + } + } // iRepeat loop + } +}// end namespace diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/InDetDDAthenaComps.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/InDetDDAthenaComps.cxx old mode 100644 new mode 100755 index 0895d8416b747ee16df2f227717e9b4a95dbc679..a85d6e394683d47b873507ba56869549d4577329 --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/InDetDDAthenaComps.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/InDetDDAthenaComps.cxx @@ -1,41 +1,35 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/InDetDDAthenaComps.h" namespace InDetDD { - -AthenaComps::AthenaComps(const std::string & msgStreamName) - : m_msg(msgStreamName), + AthenaComps::AthenaComps(const std::string& msgStreamName) + : m_msg(msgStreamName), m_detStore(0), m_geoDbTagSvc(0), m_rdbAccessSvc(0), m_geometryDBSvc(0) -{} - -void -AthenaComps::setDetStore(StoreGateSvc * detStore) -{ - m_detStore = detStore; -} - -void -AthenaComps::setGeoDbTagSvc(IGeoDbTagSvc * geoDbTagSvc) -{ - m_geoDbTagSvc = geoDbTagSvc; -} - -void -AthenaComps::setRDBAccessSvc(IRDBAccessSvc * rdbAccessSvc) -{ - m_rdbAccessSvc = rdbAccessSvc; -} - -void -AthenaComps::setGeometryDBSvc(IGeometryDBSvc * geometryDBSvc) -{ - m_geometryDBSvc = geometryDBSvc; -} - + {} + + void + AthenaComps::setDetStore(StoreGateSvc* detStore) { + m_detStore = detStore; + } + + void + AthenaComps::setGeoDbTagSvc(IGeoDbTagSvc* geoDbTagSvc) { + m_geoDbTagSvc = geoDbTagSvc; + } + + void + AthenaComps::setRDBAccessSvc(IRDBAccessSvc* rdbAccessSvc) { + m_rdbAccessSvc = rdbAccessSvc; + } + + void + AthenaComps::setGeometryDBSvc(IGeometryDBSvc* geometryDBSvc) { + m_geometryDBSvc = geometryDBSvc; + } } diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/InDetMaterialManager.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/InDetMaterialManager.cxx index 4379f3c34d789ebb599c9f6fd4a428fa9cf8b8f8..f8ae66df9ce4848ed96b55a164611b6d00a5d8c1 100755 --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/InDetMaterialManager.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/InDetMaterialManager.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/InDetMaterialManager.h" #include "InDetGeoModelUtils/InDetDDAthenaComps.h" @@ -19,28 +19,26 @@ #include <cmath> #include <stdexcept> -// Constructor -InDetMaterialManager::InDetMaterialManager(const std::string & managerName, - StoreGateSvc* detStore) +// Constructor +InDetMaterialManager::InDetMaterialManager(const std::string& managerName, + StoreGateSvc* detStore) : m_managerName(managerName), - m_msg(managerName), - m_extraFunctionality(false), - m_athenaComps(0) -{ + m_msg(managerName), + m_extraFunctionality(false), + m_athenaComps(0) { m_materialManager = retrieveManager(detStore); -} - -// Constructor -InDetMaterialManager::InDetMaterialManager(const std::string & managerName, - StoreGateSvc* detStore, - IRDBRecordset_ptr weightTable, - const std::string & space, - bool extraFunctionality) +} + +// Constructor +InDetMaterialManager::InDetMaterialManager(const std::string& managerName, + StoreGateSvc* detStore, + IRDBRecordset_ptr weightTable, + const std::string& space, + bool extraFunctionality) : m_managerName(managerName), - m_msg(managerName), - m_extraFunctionality(extraFunctionality), - m_athenaComps(0) -{ + m_msg(managerName), + m_extraFunctionality(extraFunctionality), + m_athenaComps(0) { m_materialManager = retrieveManager(detStore); if (weightTable) addWeightTable(weightTable, space); @@ -50,94 +48,81 @@ InDetMaterialManager::InDetMaterialManager(const std::string & managerName, //m_weightMap["sct::CoolingBlock"] = MaterialByWeight(2*CLHEP::gram); //m_weightMap["sct::BrlHybrid"] = MaterialByWeight("sct::Hybrid", 8*CLHEP::gram); //m_weightMap["sct::FwdHybrid"] = MaterialByWeight("std::Carbon", 7.662*CLHEP::gram); +} -} - -InDetMaterialManager::InDetMaterialManager(const std::string & managerName, StoreGateSvc* detStore, - IRDBRecordset_ptr weightTable, - IRDBRecordset_ptr compositionTable, - const std::string & space) +InDetMaterialManager::InDetMaterialManager(const std::string& managerName, StoreGateSvc* detStore, + IRDBRecordset_ptr weightTable, + IRDBRecordset_ptr compositionTable, + const std::string& space) : m_managerName(managerName), - m_msg(managerName), - m_extraFunctionality(true), - m_athenaComps(0) -{ - m_materialManager = retrieveManager(detStore); - - if (weightTable) addWeightTable(weightTable, space); - if (compositionTable) addCompositionTable(compositionTable, space); -} + m_msg(managerName), + m_extraFunctionality(true), + m_athenaComps(0) { + m_materialManager = retrieveManager(detStore); + if (weightTable) addWeightTable(weightTable, space); + if (compositionTable) addCompositionTable(compositionTable, space); +} -InDetMaterialManager::InDetMaterialManager(const std::string & managerName, - const InDetDD::AthenaComps * athenaComps) +InDetMaterialManager::InDetMaterialManager(const std::string& managerName, + const InDetDD::AthenaComps* athenaComps) : m_managerName(managerName), - m_msg(managerName), - m_extraFunctionality(true), - m_athenaComps(athenaComps) -{ + m_msg(managerName), + m_extraFunctionality(true), + m_athenaComps(athenaComps) { m_materialManager = retrieveManager(athenaComps->detStore()); addTextFileMaterials(); } - -InDetMaterialManager::~InDetMaterialManager() -{ +InDetMaterialManager::~InDetMaterialManager() { // Dereference the materials. MaterialStore::const_iterator iter; for (iter = m_store.begin(); iter != m_store.end(); ++iter) { iter->second->unref(); } - } -const AbsMaterialManager * -InDetMaterialManager::retrieveManager(StoreGateSvc* detStore) -{ - - const StoredMaterialManager * theGeoMaterialManager = nullptr; - +const AbsMaterialManager* +InDetMaterialManager::retrieveManager(StoreGateSvc* detStore) { + const StoredMaterialManager* theGeoMaterialManager = nullptr; + if (StatusCode::SUCCESS != detStore->retrieve(theGeoMaterialManager, "MATERIALS")) { msg(MSG::FATAL) << "Cannot locate Materials"; } - - return theGeoMaterialManager; + return theGeoMaterialManager; } - -const GeoElement* -InDetMaterialManager::getElement(const std::string & elementName) const -{ +const GeoElement* +InDetMaterialManager::getElement(const std::string& elementName) const { return m_materialManager->getElement(elementName); } -const GeoMaterial* -InDetMaterialManager::getMaterial(const std::string & materialName) -{ +const GeoMaterial* +InDetMaterialManager::getMaterial(const std::string& materialName) { return extraScaledMaterial(materialName, getMaterialInternal(materialName)); } -bool InDetMaterialManager::hasMaterial(const std::string &materialName) const { - return m_store.find(materialName) != m_store.end(); +bool +InDetMaterialManager::hasMaterial(const std::string& materialName) const { + return m_store.find(materialName) != m_store.end(); } -const GeoMaterial* -InDetMaterialManager::getMaterialInternal(const std::string & materialName) const -{ - // First check local store of materials. If not found then get it from the GeoModel +const GeoMaterial* +InDetMaterialManager::getMaterialInternal(const std::string& materialName) const { + // First check local store of materials. If not found then get it from the GeoModel // manager. const GeoMaterial* material = getAdditionalMaterial(materialName); + if (!material) { // This prints error message if not found. material = m_materialManager->getMaterial(materialName); - } + } return material; } -const GeoMaterial* -InDetMaterialManager::getAdditionalMaterial(const std::string & materialName) const -{ +const GeoMaterial* +InDetMaterialManager::getAdditionalMaterial(const std::string& materialName) const { MaterialStore::const_iterator iter; if ((iter = m_store.find(materialName)) != m_store.end()) { return iter->second; @@ -146,71 +131,67 @@ InDetMaterialManager::getAdditionalMaterial(const std::string & materialName) co } } -const GeoMaterial * -InDetMaterialManager::getCompositeMaterialForVolume(const std::string & newMatName, - const double volumeTot, - const double volume1, const std::string & matName1, - const double volume2, const std::string & matName2 - ) -{ - +const GeoMaterial* +InDetMaterialManager::getCompositeMaterialForVolume(const std::string& newMatName, + const double volumeTot, + const double volume1, const std::string& matName1, + const double volume2, const std::string& matName2 + ) { std::vector<std::string> baseMaterials; std::vector<double> fracWeight; baseMaterials.reserve(2); fracWeight.reserve(2); - msg(MSG::DEBUG)<<"Composite material : "<<volumeTot/CLHEP::cm3<<" = "<<volume1/CLHEP::cm3<<" + "<<volume2/CLHEP::cm3<<endmsg; - msg(MSG::DEBUG)<<"Composite material : "<<matName1<<" "<<matName2<<endmsg; - + msg(MSG::DEBUG) << "Composite material : " << volumeTot / CLHEP::cm3 << " = " << volume1 / CLHEP::cm3 << " + " << + volume2 / CLHEP::cm3 << endmsg; + msg(MSG::DEBUG) << "Composite material : " << matName1 << " " << matName2 << endmsg; + double density1, density2; - + MaterialWeightMap::const_iterator iter; if ((iter = m_weightMap.find(matName1)) != m_weightMap.end()) { - const GeoMaterial* mat1 = getMaterialForVolume(matName1,volume1); - density1=mat1->getDensity(); - msg(MSG::DEBUG)<<"Composite material 1 - weight : "<<density1/(CLHEP::gram/CLHEP::cm3)<<endmsg; - } - else { + const GeoMaterial* mat1 = getMaterialForVolume(matName1, volume1); + density1 = mat1->getDensity(); + msg(MSG::DEBUG) << "Composite material 1 - weight : " << density1 / (CLHEP::gram / CLHEP::cm3) << endmsg; + } else { const GeoMaterial* mat1 = getMaterial(matName1); - density1=mat1->getDensity(); - msg(MSG::DEBUG)<<"Composite material 1 - standard : "<<density1/(CLHEP::gram/CLHEP::cm3)<<endmsg; - } + density1 = mat1->getDensity(); + msg(MSG::DEBUG) << "Composite material 1 - standard : " << density1 / (CLHEP::gram / CLHEP::cm3) << endmsg; + } if ((iter = m_weightMap.find(matName2)) != m_weightMap.end()) { - const GeoMaterial* mat2 = getMaterialForVolume(matName2,volume2); - density2=mat2->getDensity(); - msg(MSG::DEBUG)<<"Composite material 2 - weight : "<<density2/(CLHEP::gram/CLHEP::cm3)<<endmsg; - } - else { + const GeoMaterial* mat2 = getMaterialForVolume(matName2, volume2); + density2 = mat2->getDensity(); + msg(MSG::DEBUG) << "Composite material 2 - weight : " << density2 / (CLHEP::gram / CLHEP::cm3) << endmsg; + } else { const GeoMaterial* mat2 = getMaterial(matName2); - density2=mat2->getDensity(); - msg(MSG::DEBUG)<<"Composite material 2 - standard : "<<density2/(CLHEP::gram/CLHEP::cm3)<<endmsg; - } + density2 = mat2->getDensity(); + msg(MSG::DEBUG) << "Composite material 2 - standard : " << density2 / (CLHEP::gram / CLHEP::cm3) << endmsg; + } - double weight1=density1*volume1; - double weight2=density2*volume2; - double invWeightTot=1.0/(weight1+weight2); + double weight1 = density1 * volume1; + double weight2 = density2 * volume2; + double invWeightTot = 1.0 / (weight1 + weight2); - double density=(weight1+weight2)/volumeTot; + double density = (weight1 + weight2) / volumeTot; - double frac1=weight1/(weight1+weight2); - double frac2=weight2/(weight1+weight2); - double density_2=1.0/(frac1/density1+frac2/density2); - double density_3=(weight1+weight2)/(volume1+volume2); - msg(MSG::DEBUG)<<"-> weights : "<<weight1/(CLHEP::gram)<<" "<<weight2/(CLHEP::gram)<<endmsg; - msg(MSG::DEBUG)<<"-> density : "<<density/(CLHEP::gram/CLHEP::cm3)<<" "<<density_2/(CLHEP::gram/CLHEP::cm3)<<" "<<density_3/(CLHEP::gram/CLHEP::cm3)<<endmsg; + double frac1 = weight1 / (weight1 + weight2); + double frac2 = weight2 / (weight1 + weight2); + double density_2 = 1.0 / (frac1 / density1 + frac2 / density2); + double density_3 = (weight1 + weight2) / (volume1 + volume2); + msg(MSG::DEBUG) << "-> weights : " << weight1 / (CLHEP::gram) << " " << weight2 / (CLHEP::gram) << endmsg; + msg(MSG::DEBUG) << "-> density : " << density / (CLHEP::gram / CLHEP::cm3) << " " << density_2 / + (CLHEP::gram / CLHEP::cm3) << " " << density_3 / (CLHEP::gram / CLHEP::cm3) << endmsg; baseMaterials.push_back(matName1); baseMaterials.push_back(matName2); - fracWeight.push_back(weight1*invWeightTot); - fracWeight.push_back(weight2*invWeightTot); + fracWeight.push_back(weight1 * invWeightTot); + fracWeight.push_back(weight2 * invWeightTot); - return getMaterial(newMatName,baseMaterials,fracWeight,density); + return getMaterial(newMatName, baseMaterials, fracWeight, density); } - - // This creates a new material with specified density. // If a newName is supplied it creates the new material even if the orginal material @@ -222,107 +203,101 @@ InDetMaterialManager::getCompositeMaterialForVolume(const std::string & newMatNa // name. -const GeoMaterial* -InDetMaterialManager::getMaterial(const std::string & origMaterialName, - double density, - const std::string & newName) -{ - return extraScaledMaterial(origMaterialName, newName, - getMaterialInternal(origMaterialName, density, newName)); +const GeoMaterial* +InDetMaterialManager::getMaterial(const std::string& origMaterialName, + double density, + const std::string& newName) { + return extraScaledMaterial(origMaterialName, newName, + getMaterialInternal(origMaterialName, density, newName)); } -const GeoMaterial* -InDetMaterialManager::getMaterialInternal(const std::string & origMaterialName, - double density, - const std::string & newName) -{ +const GeoMaterial* +InDetMaterialManager::getMaterialInternal(const std::string& origMaterialName, + double density, + const std::string& newName) { std::string newName2 = newName; bool newNameProvided = !newName2.empty(); if (!newNameProvided) { - newName2 = origMaterialName+"Modified"; + newName2 = origMaterialName + "Modified"; } - const GeoMaterial * newMaterial = 0; - + const GeoMaterial* newMaterial = 0; + // First see if we already have the modified material const GeoMaterial* material = getAdditionalMaterial(newName2); if (material) { if (!compareDensity(material->getDensity(), density)) { - msg(MSG::WARNING) << "Density is not consistent for material " << newName2 - << " "<<material->getDensity()/(CLHEP::gram/CLHEP::cm3) - <<" / "<<density/(CLHEP::gram/CLHEP::cm3)<<endmsg; - } + msg(MSG::WARNING) << "Density is not consistent for material " << newName2 + << " " << material->getDensity() / (CLHEP::gram / CLHEP::cm3) + << " / " << density / (CLHEP::gram / CLHEP::cm3) << endmsg; + } newMaterial = material; } else { - const GeoMaterial * origMaterial = getMaterialInternal(origMaterialName); + const GeoMaterial* origMaterial = getMaterialInternal(origMaterialName); newMaterial = origMaterial; if (origMaterial) { // If no new name was provided we check if the density is compatible // and if so we return the original material. if (newNameProvided || !compareDensity(origMaterial->getDensity(), density)) { - // create new material - GeoMaterial * newMaterialTmp = new GeoMaterial(newName2, density); - newMaterialTmp->add(const_cast<GeoMaterial*>(origMaterial), 1.); - addMaterial(newMaterialTmp); - newMaterial = newMaterialTmp; + // create new material + GeoMaterial* newMaterialTmp = new GeoMaterial(newName2, density); + newMaterialTmp->add(const_cast<GeoMaterial*>(origMaterial), 1.); + addMaterial(newMaterialTmp); + newMaterial = newMaterialTmp; } } } return newMaterial; } -const GeoMaterial* -InDetMaterialManager::getMaterialScaled(const std::string & origMaterialName, - double scaleFactor, - const std::string & newName) -{ - return extraScaledMaterial(origMaterialName, newName, - getMaterialScaledInternal(origMaterialName, scaleFactor, newName)); +const GeoMaterial* +InDetMaterialManager::getMaterialScaled(const std::string& origMaterialName, + double scaleFactor, + const std::string& newName) { + return extraScaledMaterial(origMaterialName, newName, + getMaterialScaledInternal(origMaterialName, scaleFactor, newName)); } -const GeoMaterial* -InDetMaterialManager::getMaterialScaledInternal(const std::string & origMaterialName, - double scaleFactor, - const std::string & newName) -{ +const GeoMaterial* +InDetMaterialManager::getMaterialScaledInternal(const std::string& origMaterialName, + double scaleFactor, + const std::string& newName) { // Don't allow large scale factors if (scaleFactor > 1000 || scaleFactor < 0.001) { msg(MSG::ERROR) << "Scale factor must be between 0.001 and 1000." << endmsg; return 0; } - const GeoMaterial * origMaterial = getMaterialInternal(origMaterialName); + const GeoMaterial* origMaterial = getMaterialInternal(origMaterialName); - // If scalefactor is 1 and no new name is requested + // If scalefactor is 1 and no new name is requested // then just return the orginal material if (newName.empty() && scaleFactor == 1.) return origMaterial; - const GeoMaterial * newMaterial = 0; - + const GeoMaterial* newMaterial = 0; + if (origMaterial) { double density = origMaterial->getDensity() * scaleFactor; std::string newName2 = newName; if (newName2.empty()) { // Create name using the scale factor. - int scaleInt = static_cast<int>(scaleFactor*10000); - int scale1 = scaleInt/10000; - int scale2 = scaleInt%10000; + int scaleInt = static_cast<int>(scaleFactor * 10000); + int scale1 = scaleInt / 10000; + int scale2 = scaleInt % 10000; std::ostringstream os; os << origMaterialName << scale1 << "_" << std::setw(4) << std::setfill('0') << scale2; newName2 = os.str(); } - + newMaterial = getMaterialInternal(origMaterialName, density, newName2); } return newMaterial; } - -void -InDetMaterialManager::addMaterial(GeoMaterial* material) -{ +void +InDetMaterialManager::addMaterial(GeoMaterial* material) { std::string name(material->getName()); if (m_store.find(name) != m_store.end()) { msg(MSG::WARNING) << "Ignoring attempt to redefine an existing material: " << name << endmsg; @@ -331,25 +306,22 @@ InDetMaterialManager::addMaterial(GeoMaterial* material) material->unref(); //std::cout << m_store[name] << std::endl; } else { - material->lock(); + material->lock(); material->ref(); m_store[name] = material; - - if (msgLvl(MSG::DEBUG)) - msg(MSG::DEBUG) << "Created new material: " << name << ", " << material->getDensity()/(CLHEP::g/CLHEP::cm3) << " CLHEP::g/CLHEP::cm3" << endmsg; + if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Created new material: " << name << ", " << material->getDensity() / + (CLHEP::g / CLHEP::cm3) << " CLHEP::g/CLHEP::cm3" << endmsg; } } -bool -InDetMaterialManager::compareDensity(double d1, double d2) const -{ - return (std::abs(d1/d2 - 1.) < 1e-5); +bool +InDetMaterialManager::compareDensity(double d1, double d2) const { + return(std::abs(d1 / d2 - 1.) < 1e-5); } void -InDetMaterialManager::addWeightTable(IRDBRecordset_ptr weightTable, const std::string & space) -{ +InDetMaterialManager::addWeightTable(IRDBRecordset_ptr weightTable, const std::string& space) { if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Reading in weight table: " << weightTable->nodeName() << endmsg; // If not using geometryDBSvc revert to old version if (!db()) { @@ -358,62 +330,60 @@ InDetMaterialManager::addWeightTable(IRDBRecordset_ptr weightTable, const std::s return; } for (unsigned int i = 0; i < db()->getTableSize(weightTable); i++) { - std::string materialName = db()->getString(weightTable,"MATERIAL",i); + std::string materialName = db()->getString(weightTable, "MATERIAL", i); if (!space.empty()) { - materialName = space+"::"+materialName; + materialName = space + "::" + materialName; } std::string materialBase; - if (db()->testField(weightTable,"BASEMATERIAL",i)) { - materialBase = db()->getString(weightTable,"BASEMATERIAL",i); + if (db()->testField(weightTable, "BASEMATERIAL", i)) { + materialBase = db()->getString(weightTable, "BASEMATERIAL", i); } - double weight = db()->getDouble(weightTable,"WEIGHT",i) * CLHEP::gram; + double weight = db()->getDouble(weightTable, "WEIGHT", i) * CLHEP::gram; //std::cout << materialName << " " << materialBase << " " << weight/CLHEP::g << std::endl; bool linearWeightFlag = false; - if (m_extraFunctionality && db()->testField(weightTable,"LINWEIGHTFLAG",i)) { - linearWeightFlag = db()->getInt(weightTable,"LINWEIGHTFLAG",i); + if (m_extraFunctionality && db()->testField(weightTable, "LINWEIGHTFLAG", i)) { + linearWeightFlag = db()->getInt(weightTable, "LINWEIGHTFLAG", i); } if (m_weightMap.find(materialName) != m_weightMap.end()) { msg(MSG::WARNING) << "Material: " << materialName << " already exists in weight table" << endmsg; } else { - msg(MSG::DEBUG) << "Adding " << materialName - << " weight " << weight - << " linearWeightFlag " << linearWeightFlag - << " raw weight " << db()->getDouble(weightTable,"WEIGHT",i) - << " m_extraFunctionality " << m_extraFunctionality - << " to weight table" << endmsg; + msg(MSG::DEBUG) << "Adding " << materialName + << " weight " << weight + << " linearWeightFlag " << linearWeightFlag + << " raw weight " << db()->getDouble(weightTable, "WEIGHT", i) + << " m_extraFunctionality " << m_extraFunctionality + << " to weight table" << endmsg; m_weightMap[materialName] = MaterialByWeight(materialBase, weight, linearWeightFlag); } } } - void -InDetMaterialManager::addWeightMaterial(std::string materialName, std::string materialBase, double weight, int linearWeightFlag) -{ +InDetMaterialManager::addWeightMaterial(std::string materialName, std::string materialBase, double weight, + int linearWeightFlag) { // Weight in gr weight = weight * CLHEP::gram; if (m_weightMap.find(materialName) != m_weightMap.end()) { msg(MSG::WARNING) << "Material: " << materialName << " already exists in weight table" << endmsg; } else { - msg(MSG::DEBUG) << "Adding " << materialName - << " weight " << weight - << " linearWeightFlag " << linearWeightFlag - << " to weight table" << endmsg; + msg(MSG::DEBUG) << "Adding " << materialName + << " weight " << weight + << " linearWeightFlag " << linearWeightFlag + << " to weight table" << endmsg; m_weightMap[materialName] = MaterialByWeight(materialBase, weight, linearWeightFlag); } } void -InDetMaterialManager::addWeightTableOld(IRDBRecordset_ptr weightTable, const std::string & space) -{ +InDetMaterialManager::addWeightTableOld(IRDBRecordset_ptr weightTable, const std::string& space) { for (unsigned int i = 0; i < weightTable->size(); i++) { - const IRDBRecord * record = (*weightTable)[i]; - std::string materialName = record->getString("MATERIAL"); + const IRDBRecord* record = (*weightTable)[i]; + std::string materialName = record->getString("MATERIAL"); if (!space.empty()) { - materialName = space+"::"+materialName; + materialName = space + "::" + materialName; } std::string materialBase; if (!record->isFieldNull("BASEMATERIAL")) { @@ -436,60 +406,63 @@ InDetMaterialManager::addWeightTableOld(IRDBRecordset_ptr weightTable, const std } void -InDetMaterialManager::addCompositionTable(IRDBRecordset_ptr compositionTable, const std::string & space) -{ +InDetMaterialManager::addCompositionTable(IRDBRecordset_ptr compositionTable, const std::string& space) { if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Reading in composition table: " << compositionTable->nodeName() << endmsg; + if (!db()) { msg(MSG::ERROR) << "GeometryDBSvc not available. Unable to read in composition table." << endmsg; } for (unsigned int i = 0; i < db()->getTableSize(compositionTable); i++) { - std::string materialName = db()->getString(compositionTable,"MATERIAL",i); + std::string materialName = db()->getString(compositionTable, "MATERIAL", i); if (!space.empty()) { - materialName = space+"::"+materialName; + materialName = space + "::" + materialName; } - - std::string componentName = db()->getString(compositionTable,"COMPONENT",i); - int count = db()->getInt(compositionTable,"COUNT",i); - double factor = db()->getDouble(compositionTable,"FACTOR",i); - double actualLength = db()->getDouble(compositionTable,"ACTUALLENGTH",i); - m_matCompositionMap.insert(std::pair<std::string,MaterialComponent>(materialName, MaterialComponent(componentName, count*factor, actualLength))); + std::string componentName = db()->getString(compositionTable, "COMPONENT", i); + int count = db()->getInt(compositionTable, "COUNT", i); + double factor = db()->getDouble(compositionTable, "FACTOR", i); + double actualLength = db()->getDouble(compositionTable, "ACTUALLENGTH", i); + + m_matCompositionMap.insert(std::pair<std::string, MaterialComponent>(materialName, + MaterialComponent(componentName, + count * factor, + actualLength))); } } void -InDetMaterialManager::addScalingTable(IRDBRecordset_ptr scalingTable) -{ +InDetMaterialManager::addScalingTable(IRDBRecordset_ptr scalingTable) { if (!scalingTable) return; + if (db()->getTableSize(scalingTable) == 0) return; - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Reading in extra material scaling table: " << scalingTable->nodeName() << endmsg; + if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Reading in extra material scaling table: " << scalingTable->nodeName() << + endmsg; if (!db()) { msg(MSG::ERROR) << "GeometryDBSvc not available. Unable to read in scaling table." << endmsg; } for (unsigned int i = 0; i < db()->getTableSize(scalingTable); i++) { - std::string materialName = db()->getString(scalingTable,"MATERIAL",i); - double scalingFactor = db()->getDouble(scalingTable,"FACTOR",i); + std::string materialName = db()->getString(scalingTable, "MATERIAL", i); + double scalingFactor = db()->getDouble(scalingTable, "FACTOR", i); if (msgLvl(MSG::DEBUG)) { if (scalingFactor >= 0 || scalingFactor == 1) { - msg(MSG::DEBUG) << "Material " << materialName << " will be scaled by: " << scalingFactor << endmsg; + msg(MSG::DEBUG) << "Material " << materialName << " will be scaled by: " << scalingFactor << endmsg; } else { - // -ve or scalefactor = 1 means will not be scaled. - msg(MSG::DEBUG) << "Material " << materialName << " will be NOT be scaled." << endmsg; - } + // -ve or scalefactor = 1 means will not be scaled. + msg(MSG::DEBUG) << "Material " << materialName << " will be NOT be scaled." << endmsg; + } } if (m_scalingMap.find(materialName) != m_scalingMap.end()) { - msg(MSG::WARNING) << "Overriding material: " << materialName << " which already exists in scaling table" << endmsg; - } + msg(MSG::WARNING) << "Overriding material: " << materialName << " which already exists in scaling table" << + endmsg; + } m_scalingMap[materialName] = scalingFactor; } } - -const GeoMaterial * -InDetMaterialManager::getMaterialForVolume(const std::string & materialName, double volume, const std::string & newName) -{ +const GeoMaterial* +InDetMaterialManager::getMaterialForVolume(const std::string& materialName, double volume, const std::string& newName) { // Make sure we have a valid volume size. if (volume <= 0) { msg(MSG::ERROR) << "Invalid volume : " << volume << endmsg; @@ -508,50 +481,52 @@ InDetMaterialManager::getMaterialForVolume(const std::string & materialName, dou MaterialWeightMap::const_iterator iter; if ((iter = m_weightMap.find(materialName)) != m_weightMap.end()) { - const std::string & materialBase = iter->second.name; + const std::string& materialBase = iter->second.name; double weight = iter->second.weight; - double density = weight/volume; + double density = weight / volume; if (iter->second.linearWeightFlag) { - msg(MSG::ERROR) << "Material defined by linear weight cannot be created with getMaterialForVolume method: " << materialName << endmsg; - } - - if (msgLvl(MSG::VERBOSE)) { - msg(MSG::VERBOSE) - << "Found material in weight table - name, base, weight(CLHEP::g), volume(CLHEP::cm3), density(CLHEP::g/CLHEP::cm3): " - << materialName << ", " - << materialBase << ", " - << weight/CLHEP::gram << ", " - << volume/CLHEP::cm3 << ", " - << density/(CLHEP::g/CLHEP::cm3) << endmsg; + msg(MSG::ERROR) << "Material defined by linear weight cannot be created with getMaterialForVolume method: " << + materialName << endmsg; + } + + if (msgLvl(MSG::VERBOSE)) { + msg(MSG::VERBOSE) + << + "Found material in weight table - name, base, weight(CLHEP::g), volume(CLHEP::cm3), density(CLHEP::g/CLHEP::cm3): " + << materialName << ", " + << materialBase << ", " + << weight / CLHEP::gram << ", " + << volume / CLHEP::cm3 << ", " + << density / (CLHEP::g / CLHEP::cm3) << endmsg; } if (materialBase.empty()) { return getMaterial(materialName, density, newName); } else { if (newName.empty()) { - return getMaterial(materialBase, density, materialName); + return getMaterial(materialBase, density, materialName); } else { - return getMaterial(materialBase, density, newName); - } + return getMaterial(materialBase, density, newName); + } } } else { // If not in the weight table we just return the material. // This is not an error. - if (msgLvl(MSG::VERBOSE)) - msg(MSG::VERBOSE) - << "Material not in weight table, using standard material: " - << materialName - << ", volume(CLHEP::cm3) = " << volume/CLHEP::cm3 - << endmsg; + if (msgLvl(MSG::VERBOSE)) + msg(MSG::VERBOSE) + << "Material not in weight table, using standard material: " + << materialName + << ", volume(CLHEP::cm3) = " << volume / CLHEP::cm3 + << endmsg; return getMaterial(materialName); - } + } } - -const GeoMaterial * -InDetMaterialManager::getMaterialForVolumeLength(const std::string & materialName, double volume, double length, const std::string & newName) -{ - // In the case there is no material composition table (MaterialCompositionMap) and no linear weights are used this will +const GeoMaterial* +InDetMaterialManager::getMaterialForVolumeLength(const std::string& materialName, double volume, double length, + const std::string& newName) { + // In the case there is no material composition table (MaterialCompositionMap) and no linear weights are used this + // will // behave the same way as getMaterialForVolume. // If the material is in the MaterialCompositionMap it will build a material using the components // from that table. If any components are defined as a linear weight the length is used to calculate the @@ -566,38 +541,37 @@ InDetMaterialManager::getMaterialForVolumeLength(const std::string & materialNam } // Make sure we have a valid volume size. - if (volume <= 0 || length <= 0) { + if (volume <= 0 || length <= 0) { msg(MSG::ERROR) << "Invalid volume or length : " << volume << ", " << length << endmsg; return 0; } // First look in the predefinded collections std::pair<MaterialCompositionMap::const_iterator, MaterialCompositionMap::const_iterator> iterRange; - iterRange = m_matCompositionMap.equal_range(materialName); - if (iterRange.first != m_matCompositionMap.end()) { - - if (msgLvl(MSG::VERBOSE)) { - msg(MSG::VERBOSE) - << "Found material in material composition table:" << materialName <<endmsg; + iterRange = m_matCompositionMap.equal_range(materialName); + if (iterRange.first != m_matCompositionMap.end()) { + if (msgLvl(MSG::VERBOSE)) { + msg(MSG::VERBOSE) + << "Found material in material composition table:" << materialName << endmsg; } std::vector<double> factors; std::vector<std::string> components; for (MaterialCompositionMap::const_iterator iter = iterRange.first; iter != iterRange.second; iter++) { double factorTmp = iter->second.factor; - if (iter->second.actualLength > 0) factorTmp *= iter->second.actualLength/length; + if (iter->second.actualLength > 0) factorTmp *= iter->second.actualLength / length; factors.push_back(factorTmp); components.push_back(iter->second.name); } - return getMaterialForVolumeLength(name,components,factors,volume,length); - } + return getMaterialForVolumeLength(name, components, factors, volume, length); + } // Next look in weight table MaterialWeightMap::const_iterator iter; if ((iter = m_weightMap.find(materialName)) != m_weightMap.end()) { - const std::string & materialBase = iter->second.name; + const std::string& materialBase = iter->second.name; double weight = iter->second.weight; - double density = weight/volume; + double density = weight / volume; if (iter->second.linearWeightFlag) weight *= length; if (materialBase.empty()) { @@ -608,47 +582,44 @@ InDetMaterialManager::getMaterialForVolumeLength(const std::string & materialNam } else { // Otherwise we just return the material. // This is not an error. - if (msgLvl(MSG::VERBOSE)) - msg(MSG::VERBOSE) - << "Material not in weight table, using standard material: " - << materialName - << ", volume(CLHEP::cm3) = " << volume/CLHEP::cm3 - << endmsg; + if (msgLvl(MSG::VERBOSE)) + msg(MSG::VERBOSE) + << "Material not in weight table, using standard material: " + << materialName + << ", volume(CLHEP::cm3) = " << volume / CLHEP::cm3 + << endmsg; return getMaterial(materialName); - } + } } - -const GeoMaterial * -InDetMaterialManager::getMaterialForVolumeLength(const std::string & name, - const std::string & materialComponent, - double factor, - double volume, - double length) -{ - std::vector<std::string> tmpMaterialComponents(1,materialComponent); - std::vector<double> tmpFactors(1,factor); +const GeoMaterial* +InDetMaterialManager::getMaterialForVolumeLength(const std::string& name, + const std::string& materialComponent, + double factor, + double volume, + double length) { + std::vector<std::string> tmpMaterialComponents(1, materialComponent); + std::vector<double> tmpFactors(1, factor); return getMaterialForVolumeLength(name, tmpMaterialComponents, tmpFactors, volume, length); } -const GeoMaterial * -InDetMaterialManager::getMaterialForVolumeLength(const std::string & name, - const std::vector<std::string> & materialComponents, - const std::vector<double> factors, - double volume, - double length) -{ - +const GeoMaterial* +InDetMaterialManager::getMaterialForVolumeLength(const std::string& name, + const std::vector<std::string>& materialComponents, + const std::vector<double> factors, + double volume, + double length) { // Make sure we have a valid volume size. - if (volume <= 0 || length <= 0) { + if (volume <= 0 || length <= 0) { msg(MSG::ERROR) << "Invalid volume or length : " << volume << ", " << length << endmsg; return 0; } if (!factors.empty() && factors.size() < materialComponents.size()) { - msg(MSG::WARNING) << "getMaterialForVolumeLength: factor vector size too small. Setting remaining factors to 1." << endmsg; + msg(MSG::WARNING) << "getMaterialForVolumeLength: factor vector size too small. Setting remaining factors to 1." << + endmsg; } - + std::vector<std::string> baseMaterials; std::vector<double> fracWeight; baseMaterials.reserve(materialComponents.size()); @@ -656,48 +627,45 @@ InDetMaterialManager::getMaterialForVolumeLength(const std::string & name, double totWeight = 0; for (unsigned int iComp = 0; iComp < materialComponents.size(); ++iComp) { - - const std::string & materialName = materialComponents[iComp]; + const std::string& materialName = materialComponents[iComp]; // First search in MaterialWeightMap MaterialWeightMap::const_iterator iter; if ((iter = m_weightMap.find(materialName)) != m_weightMap.end()) { - - const std::string & materialBase = iter->second.name; + const std::string& materialBase = iter->second.name; double weight = iter->second.weight; - if (iComp < factors.size()){ - weight *= factors[iComp]; + if (iComp < factors.size()) { + weight *= factors[iComp]; } - msg(MSG::DEBUG) << "Material " << materialName - << " found in weight table, weight " << iter->second.weight / CLHEP::gram - << " factor " << factors[iComp] - << " w*fac*len " << weight*length / CLHEP::gram - << " basMat " << materialBase - << " linear? " << iter->second.linearWeightFlag << endmsg; + msg(MSG::DEBUG) << "Material " << materialName + << " found in weight table, weight " << iter->second.weight / CLHEP::gram + << " factor " << factors[iComp] + << " w*fac*len " << weight * length / CLHEP::gram + << " basMat " << materialBase + << " linear? " << iter->second.linearWeightFlag << endmsg; if (iter->second.linearWeightFlag) weight *= length; if (materialBase.empty()) { - // If no base material then name should refer to an already defined material - baseMaterials.push_back(materialName); + // If no base material then name should refer to an already defined material + baseMaterials.push_back(materialName); } else { - baseMaterials.push_back(materialBase); + baseMaterials.push_back(materialBase); } fracWeight.push_back(weight); // Will be normalized later. totWeight += weight; - } else { // If not in the weight table we look for a regular material. // I don't think this would normally be intentional so we give a warning message. /* - if (msgLvl(MSG::WARNING)) - msg(MSG::WARNING) - << "Component material not in weight table, using standard material: " - << materialName << " with weight= " - << factors.at(iComp) * length - << endmsg; - const GeoMaterial * material = getMaterialInternal(materialName); - */ + if (msgLvl(MSG::WARNING)) + msg(MSG::WARNING) + << "Component material not in weight table, using standard material: " + << materialName << " with weight= " + << factors.at(iComp) * length + << endmsg; + const GeoMaterial * material = getMaterialInternal(materialName); + */ // In this case the factor should correspond to the linear weight double weight = factors.at(iComp) * length * CLHEP::gram; @@ -709,95 +677,93 @@ InDetMaterialManager::getMaterialForVolumeLength(const std::string & name, } } - if (msgLvl(MSG::VERBOSE)) { + if (msgLvl(MSG::VERBOSE)) { msg(MSG::VERBOSE) << "Creating material from multiple components: " << name << endmsg; for (unsigned int i = 0; i < materialComponents.size(); ++i) { - msg(MSG::VERBOSE) << " Component " << i << ": Name = " << baseMaterials[i] - << " Weight(CLHEP::g) = " << fracWeight[i]/CLHEP::g << endmsg; + msg(MSG::VERBOSE) << " Component " << i << ": Name = " << baseMaterials[i] + << " Weight(CLHEP::g) = " << fracWeight[i] / CLHEP::g << endmsg; } } - + for (unsigned int i = 0; i < fracWeight.size(); ++i) { fracWeight[i] /= totWeight; } - double density = totWeight/volume; + double density = totWeight / volume; - return getMaterial(name,baseMaterials,fracWeight,density); + return getMaterial(name, baseMaterials, fracWeight, density); } - // Add materials assuming they simply occupy the same volume. /* -const GeoMaterial* -InDetMaterialManager::getMaterial(const std::vector<const GeoMaterial *> & materialComponents, - const std::string & newName) -{ - const GeoMaterial * newMaterial = 0; - std::vector<double> fracWeight; - fracWeight.reserve(materialComponents.size()); - - for (unsigned int i = 0; i < materialComponents.size(); i++) { + const GeoMaterial* + InDetMaterialManager::getMaterial(const std::vector<const GeoMaterial *> & materialComponents, + const std::string & newName) + { + const GeoMaterial * newMaterial = 0; + std::vector<double> fracWeight; + fracWeight.reserve(materialComponents.size()); + + for (unsigned int i = 0; i < materialComponents.size(); i++) { const GeoMaterial * origMaterial = materialComponents[i]; double weight = origMaterial->getDensity(); fracWeight.push_back(weight); totWeight += weight; - } - for (unsigned int i = 0; i < fracWeight.size(); ++i) { + } + for (unsigned int i = 0; i < fracWeight.size(); ++i) { fracWeight[i] /= totWeight; - } - return getMaterial(materialComponents, fracWeight, totWeight, newName); -} + } + return getMaterial(materialComponents, fracWeight, totWeight, newName); + } -const GeoMaterial* -InDetMaterialManager::getMaterial(const std::vector<std::string> & materialComponents, - const std::string & newName) -{ - const GeoMaterial * newMaterial = 0; - - // First see if we already have the modified material - const GeoMaterial* material = getAdditionalMaterial(newName); + const GeoMaterial* + InDetMaterialManager::getMaterial(const std::vector<std::string> & materialComponents, + const std::string & newName) + { + const GeoMaterial * newMaterial = 0; + + // First see if we already have the modified material + const GeoMaterial* material = getAdditionalMaterial(newName); - for (unsigned int i = 0; i < materialComponents.size(); i++) { - const GeoMaterial * origMaterial = getMaterial(materialComponents[i]); + for (unsigned int i = 0; i < materialComponents.size(); i++) { + const GeoMaterial * origMaterial = getMaterial(materialComponents[i]); components.push_back(origMaterial); - } - return getMaterial(components, newName); -} -*/ + } + return getMaterial(components, newName); + } + */ -const GeoMaterial* -InDetMaterialManager::getMaterial(const std::string & name, - const std::vector<std::string> & materialComponents, - const std::vector<double> & fracWeight, - double density) -{ +const GeoMaterial* +InDetMaterialManager::getMaterial(const std::string& name, + const std::vector<std::string>& materialComponents, + const std::vector<double>& fracWeight, + double density) { return extraScaledMaterial(name, getMaterialInternal(name, materialComponents, fracWeight, density)); } -const GeoMaterial* -InDetMaterialManager::getMaterialInternal(const std::string & name, - const std::vector<std::string> & materialComponents, - const std::vector<double> & fracWeight, - double density) -{ - const GeoMaterial * newMaterial = 0; - +const GeoMaterial* +InDetMaterialManager::getMaterialInternal(const std::string& name, + const std::vector<std::string>& materialComponents, + const std::vector<double>& fracWeight, + double density) { + const GeoMaterial* newMaterial = 0; + // First see if we already have the material const GeoMaterial* material = getAdditionalMaterial(name); + if (material) { if (!compareDensity(material->getDensity(), density)) { msg(MSG::WARNING) << "Density is not consistent for material " << name << endmsg; - } + } newMaterial = material; } else { - GeoMaterial * newMaterialTmp = new GeoMaterial(name, density); + GeoMaterial* newMaterialTmp = new GeoMaterial(name, density); for (unsigned int i = 0; i < materialComponents.size(); i++) { - const GeoMaterial * origMaterial = getMaterialInternal(materialComponents[i]); + const GeoMaterial* origMaterial = getMaterialInternal(materialComponents[i]); if (origMaterial) { - newMaterialTmp->add(const_cast<GeoMaterial*>(origMaterial), fracWeight[i]); + newMaterialTmp->add(const_cast<GeoMaterial*>(origMaterial), fracWeight[i]); } else { - msg(MSG::ERROR) << "Material component missing " << materialComponents[i] << endmsg; + msg(MSG::ERROR) << "Material component missing " << materialComponents[i] << endmsg; } } addMaterial(newMaterialTmp); @@ -806,24 +772,25 @@ InDetMaterialManager::getMaterialInternal(const std::string & name, return newMaterial; } -const IGeometryDBSvc * InDetMaterialManager::db() -{ +const IGeometryDBSvc* +InDetMaterialManager::db() { if (m_athenaComps) return m_athenaComps->geomDB(); + return 0; } -void InDetMaterialManager::addTextFileMaterials() -{ +void +InDetMaterialManager::addTextFileMaterials() { const std::string materialTable = "ExtraMaterials"; const std::string componentsTable = "ExtraMatComponents"; // Look for tables ExtraMaterials and ExtraMatComponents. - // These are text file only tables where extra materials are desired or + // These are text file only tables where extra materials are desired or // one wants to override some database ones. - if (!db() || !db()->testField("","TableSize:"+materialTable) || !db()->getTableSize(materialTable) - || !db()->testField("","TableSize:"+componentsTable) || !db()->getTableSize(componentsTable)) return; + if (!db() || !db()->testField("", "TableSize:" + materialTable) || !db()->getTableSize(materialTable) + || !db()->testField("", "TableSize:" + componentsTable) || !db()->getTableSize(componentsTable)) return; + - msg(MSG::INFO) << "Extra materials being read in from text file." << endmsg; typedef std::map<std::string, MaterialDef> MatMap; @@ -831,22 +798,23 @@ void InDetMaterialManager::addTextFileMaterials() // read in material table for (unsigned int iMat = 0; iMat < db()->getTableSize(materialTable); iMat++) { - std::string materialName = db()->getString(materialTable,"NAME",iMat); - double density = db()->getDouble(materialTable,"DENSITY",iMat)*CLHEP::g/CLHEP::cm3; - materials[materialName] = MaterialDef(materialName,density); + std::string materialName = db()->getString(materialTable, "NAME", iMat); + double density = db()->getDouble(materialTable, "DENSITY", iMat) * CLHEP::g / CLHEP::cm3; + materials[materialName] = MaterialDef(materialName, density); } - + // read in material component table for (unsigned int iComp = 0; iComp < db()->getTableSize(componentsTable); iComp++) { - std::string materialName = db()->getString(componentsTable,"NAME",iComp); - std::string compName = db()->getString(componentsTable,"COMPNAME",iComp); - double fracWeight = db()->getDouble(componentsTable,"FRACTION",iComp); - MatMap::iterator iter = materials.find(materialName); - if (iter != materials.end()) { - iter->second.addComponent(compName,fracWeight); - } else { - msg(MSG::ERROR) << "Attemp to add material component, " << compName << ", to non-existing material: " << materialName << endmsg; - } + std::string materialName = db()->getString(componentsTable, "NAME", iComp); + std::string compName = db()->getString(componentsTable, "COMPNAME", iComp); + double fracWeight = db()->getDouble(componentsTable, "FRACTION", iComp); + MatMap::iterator iter = materials.find(materialName); + if (iter != materials.end()) { + iter->second.addComponent(compName, fracWeight); + } else { + msg(MSG::ERROR) << "Attemp to add material component, " << compName << ", to non-existing material: " << + materialName << endmsg; + } } //Now create the materials @@ -859,47 +827,46 @@ void InDetMaterialManager::addTextFileMaterials() while (someUndefined && matCount != matCountLast) { matCountLast = matCount; someUndefined = false; - for (MatMap::iterator iter = materials.begin(); iter != materials.end(); ++iter) { - MaterialDef & tmpMat = iter->second; + for (MatMap::iterator iter = materials.begin(); iter != materials.end(); ++iter) { + MaterialDef& tmpMat = iter->second; if (!tmpMat.isCreated()) { - // Check if any components are materials in this table and if they are defined. - // If not flag that there are undefined materials and go to next material - bool compsDefined = true; - for (unsigned int iComp = 0; iComp < tmpMat.numComponents(); ++iComp) { - std::string compName = tmpMat.compName(iComp); - MatMap::iterator iter2 = materials.find(compName); - if (iter2 != materials.end()) { - if (!iter2->second.isCreated()) { - compsDefined = false; - break; - } - } - } - if (compsDefined) { - createMaterial(tmpMat); - tmpMat.setCreated(); - matCount++; - } else { - someUndefined = true; - } + // Check if any components are materials in this table and if they are defined. + // If not flag that there are undefined materials and go to next material + bool compsDefined = true; + for (unsigned int iComp = 0; iComp < tmpMat.numComponents(); ++iComp) { + std::string compName = tmpMat.compName(iComp); + MatMap::iterator iter2 = materials.find(compName); + if (iter2 != materials.end()) { + if (!iter2->second.isCreated()) { + compsDefined = false; + break; + } + } + } + if (compsDefined) { + createMaterial(tmpMat); + tmpMat.setCreated(); + matCount++; + } else { + someUndefined = true; + } } } } - - + + if (someUndefined) { msg(MSG::ERROR) << "Not all materials could be defined due to cyclic definitions" << endmsg; } } - -void InDetMaterialManager::createMaterial(const MaterialDef & material) -{ +void +InDetMaterialManager::createMaterial(const MaterialDef& material) { if (material.numComponents() == 0) { msg(MSG::ERROR) << "Material has no components: " << material.name() << endmsg; return; } - + // If total of fractions is greater than 1.1 then assume material is define by ratio of atoms. double totWeight = material.totalFraction(); bool byAtomicRatio = false; @@ -907,58 +874,63 @@ void InDetMaterialManager::createMaterial(const MaterialDef & material) byAtomicRatio = true; for (unsigned int i = 0; i < material.numComponents(); i++) { if (material.compName(i).find("::") != std::string::npos) { - // If component name has "::" in it then its not an element. - msg(MSG::ERROR) << "Material, " << material.name() - << ", is assumed to be defined by atomic ratio (due to total fraction > 1) but component is not an element: " - << material.compName(i) << endmsg; - return; + // If component name has "::" in it then its not an element. + msg(MSG::ERROR) << "Material, " << material.name() + << + ", is assumed to be defined by atomic ratio (due to total fraction > 1) but component is not an element: " + << material.compName(i) << endmsg; + return; } - const GeoElement * element = getElement(material.compName(i)); + const GeoElement* element = getElement(material.compName(i)); if (!element) { - msg(MSG::ERROR) << "Error making material " << material.name() << ". Element not found: " << material.compName(i) << endmsg; - return; + msg(MSG::ERROR) << "Error making material " << material.name() << ". Element not found: " << + material.compName(i) << endmsg; + return; } totWeight += material.fraction(i) * element->getA(); } } else { // Check if total fraction is close to 1. if (std::abs(totWeight - 1) > 0.01) { - msg(MSG::WARNING) << "Total fractional weight does not sum to 1. Will renormalize. Total = " << totWeight << endmsg; + msg(MSG::WARNING) << "Total fractional weight does not sum to 1. Will renormalize. Total = " << totWeight << + endmsg; } - } + } // Now build the material - GeoMaterial * newMaterial = new GeoMaterial(material.name(),material.density()); - if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Creating material: " << material.name() - << " with density: " << material.density()/(CLHEP::g/CLHEP::cm3) << endmsg; + GeoMaterial* newMaterial = new GeoMaterial(material.name(), material.density()); + if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << "Creating material: " << material.name() + << " with density: " << material.density() / (CLHEP::g / CLHEP::cm3) << + endmsg; for (unsigned int i = 0; i < material.numComponents(); i++) { double fracWeight = material.fraction(i) / totWeight; if (material.compName(i).find("::") == std::string::npos) { - const GeoElement * element = getElement(material.compName(i)); + const GeoElement* element = getElement(material.compName(i)); if (!element) { - msg(MSG::ERROR) << "Error making material " << material.name() << ". Element not found: " << material.compName(i) << endmsg; - // delete the partially created material - newMaterial->ref(); - newMaterial->unref(); - return; + msg(MSG::ERROR) << "Error making material " << material.name() << ". Element not found: " << + material.compName(i) << endmsg; + // delete the partially created material + newMaterial->ref(); + newMaterial->unref(); + return; } if (byAtomicRatio) { - fracWeight = material.fraction(i) * element->getA() / totWeight; + fracWeight = material.fraction(i) * element->getA() / totWeight; } newMaterial->add(const_cast<GeoElement*>(element), fracWeight); if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " Component: " << material.compName(i) << " " << fracWeight << endmsg; - } else { - const GeoMaterial * materialTmp = getMaterialInternal(material.compName(i)); + const GeoMaterial* materialTmp = getMaterialInternal(material.compName(i)); if (!materialTmp) { - msg(MSG::ERROR) << "Error making material " << material.name() << ". Component not found: " << material.compName(i) << endmsg; - // delete the partially created material - newMaterial->ref(); - newMaterial->unref(); - return; + msg(MSG::ERROR) << "Error making material " << material.name() << ". Component not found: " << + material.compName(i) << endmsg; + // delete the partially created material + newMaterial->ref(); + newMaterial->unref(); + return; } if (byAtomicRatio) { - // Should not happen as already checked that all components were elements. - msg(MSG::ERROR) << "Unexpected Error" << endmsg; + // Should not happen as already checked that all components were elements. + msg(MSG::ERROR) << "Unexpected Error" << endmsg; } newMaterial->add(const_cast<GeoMaterial*>(materialTmp), fracWeight); if (msgLvl(MSG::DEBUG)) msg(MSG::DEBUG) << " Component: " << material.compName(i) << " " << fracWeight << endmsg; @@ -968,35 +940,33 @@ void InDetMaterialManager::createMaterial(const MaterialDef & material) addMaterial(newMaterial); } - - InDetMaterialManager::MaterialDef::MaterialDef() : m_density(0), - m_created(false) + m_created(false) {} - -InDetMaterialManager::MaterialDef::MaterialDef(const std::string & name, double density) - : m_name(name), - m_density(density), - m_created(false) + +InDetMaterialManager::MaterialDef::MaterialDef(const std::string& name, double density) + : m_name(name), + m_density(density), + m_created(false) {} -void InDetMaterialManager::MaterialDef::addComponent(const std::string & compName, double fraction) -{ +void +InDetMaterialManager::MaterialDef::addComponent(const std::string& compName, double fraction) { m_components.push_back(compName); m_fractions.push_back(fraction); } -double InDetMaterialManager::MaterialDef::totalFraction() const -{ +double +InDetMaterialManager::MaterialDef::totalFraction() const { double sum = 0; + for (unsigned int i = 0; i < m_fractions.size(); i++) { sum += m_fractions[i]; } return sum; } - // We need the original name as the GeoMaterial from the standard // material manager has its namespace dropped. We have two versions // of extraScaledMaterial. One where two names are provided. In this @@ -1004,42 +974,41 @@ double InDetMaterialManager::MaterialDef::totalFraction() const // materialName is used. The other just has one name and that is the // one that is used. -const GeoMaterial* -InDetMaterialManager::extraScaledMaterial(const std::string & materialName, - const std::string & newName, - const GeoMaterial * origMaterial) -{ +const GeoMaterial* +InDetMaterialManager::extraScaledMaterial(const std::string& materialName, + const std::string& newName, + const GeoMaterial* origMaterial) { if (newName.empty()) { return extraScaledMaterial(materialName, origMaterial); } else { - return extraScaledMaterial(newName, origMaterial); + return extraScaledMaterial(newName, origMaterial); } } -const GeoMaterial* -InDetMaterialManager::extraScaledMaterial(const std::string & materialName, const GeoMaterial * origMaterial) -{ - if (!origMaterial) throw std::runtime_error(std::string("Invalid material: ")+materialName); +const GeoMaterial* +InDetMaterialManager::extraScaledMaterial(const std::string& materialName, const GeoMaterial* origMaterial) { + if (!origMaterial) throw std::runtime_error(std::string("Invalid material: ") + materialName); double scaleFactor = getExtraScaleFactor(materialName); // -1 (or any -ve number) indicates material is not scaled. And if the scale factor // is 1 then there is no need to create a new material. - if (scaleFactor < 0 || scaleFactor == 1 || materialName.find("Ether")!=std::string::npos) return origMaterial; + if (scaleFactor < 0 || scaleFactor == 1 || materialName.find("Ether") != std::string::npos) return origMaterial; + if (scaleFactor == 0) return getMaterialInternal("std::Vacuum"); - - std::string newName = materialName+"_ExtraScaling"; - + + std::string newName = materialName + "_ExtraScaling"; + // Check if it is already made. - const GeoMaterial * newMaterial = getAdditionalMaterial(newName); + const GeoMaterial* newMaterial = getAdditionalMaterial(newName); // Already made so we return it. if (newMaterial) return newMaterial; - - // Otherwise we need to make it. + + // Otherwise we need to make it. double density = origMaterial->getDensity() * scaleFactor; - // create new material - GeoMaterial * newMaterialTmp = new GeoMaterial(newName, density); + // create new material + GeoMaterial* newMaterialTmp = new GeoMaterial(newName, density); newMaterialTmp->add(const_cast<GeoMaterial*>(origMaterial), 1.); addMaterial(newMaterialTmp); newMaterial = newMaterialTmp; @@ -1048,26 +1017,25 @@ InDetMaterialManager::extraScaledMaterial(const std::string & materialName, cons } double -InDetMaterialManager::getExtraScaleFactor(const std::string & materialName) -{ - // If name is found in map we return the corresponding scale factor. +InDetMaterialManager::getExtraScaleFactor(const std::string& materialName) { + // If name is found in map we return the corresponding scale factor. // The special name "ALL" indicates all materials are scaled. - // Individual materials can be excluded from scaling by giving either + // Individual materials can be excluded from scaling by giving either // a -ve scaling factor or just specifying a scaling factor of 1. // A scaling factor of 0 means the material will be replaced by vacuum. ExtraScaleFactorMap::const_iterator iter = m_scalingMap.find(materialName); - if (iter != m_scalingMap.end()) { + if (iter != m_scalingMap.end()) { return iter->second; } else { // Check for special names // ALL means everything scaled. Do not scale air or vacuum (unless explicity requested) iter = m_scalingMap.find("ALL"); - if (iter != m_scalingMap.end() && materialName != "std::Air" && materialName != "std::Vacuum") { - return iter->second; + if (iter != m_scalingMap.end() && materialName != "std::Air" && materialName != "std::Vacuum") { + return iter->second; } } - + // If not found then return -1 to indicate material is not to be scaled. return -1; } diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/PairIndexMap.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/PairIndexMap.cxx old mode 100644 new mode 100755 index 4cf88367201adf13fffbda4a1a628f142e64935a..f751382680012ef5ff1b648bc5a68c3e47de4bde --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/PairIndexMap.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/PairIndexMap.cxx @@ -1,21 +1,20 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/PairIndexMap.h" namespace InDetDD { - void - PairIndexMap::add(int first, int second, int value) - { - m_map[std::make_pair(first,second)]=value; + PairIndexMap::add(int first, int second, int value) { + m_map[std::make_pair(first, second)] = value; } - - int PairIndexMap::find(int first, int second) const - { - MapType::const_iterator iter = m_map.find(std::make_pair(first,second)); + + int + PairIndexMap::find(int first, int second) const { + MapType::const_iterator iter = m_map.find(std::make_pair(first, second)); if (iter == m_map.end()) return -1; + return iter->second; } } diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ServiceVolume.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ServiceVolume.cxx index 825e9e4726a517b4ea088ba05ebda5d9c1cd95ce..43833bac89f9d8ed9881b20921e9544969ca1471 100755 --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ServiceVolume.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ServiceVolume.cxx @@ -1,25 +1,25 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ // // This class to hold general services // -// The services support several different shapes. The meaning of the parameters -// depends on the particular shape: -// -// +// The services support several different shapes. The meaning of the parameters +// depends on the particular shape: +// +// // TUBE or empty // Ignored: RIN2,ROUT2,PHI,WIDTH,REPEAT // TUBS // Ignored: RIN2,ROUT2 // PHI: phi start location of tube sector -// WIDTH (CLHEP::deg): phi width of sector +// WIDTH (CLHEP::deg): phi width of sector // REPEAT: Repeat the tube sector this many times in phi with equal distance between them. // CONS, CONE // WIDTH,REPEAT ignored if CONE // RIN2,ROUT2: rmin, rmx at zmax. Same as RIN, ROUT if <=0. -// PHI, WIDTH, REPEAT same as TUBS +// PHI, WIDTH, REPEAT same as TUBS // PGON // Ignored: WIDTH // RIN,ROUT,RIN2,ROUT2 defined at corner of polygon. @@ -52,14 +52,14 @@ // Ignored: ROUT, RIN2, ROUT2 // RIN: Radial position of center of tube // PHI: phi position of center -// WIDTH (CLHEP::mm): diameter +// WIDTH (CLHEP::mm): diameter // REPEAT: Repeat this many times in phi with equal distance between them. // ROD2 (hollow tube not centered around Z axis) // Ignored: ROUT, ROUT2 // RIN: Radial position of center of tube -// RIN2: inner radius +// RIN2: inner radius // PHI: phi position of center -// WIDTH (CLHEP::mm): diameter +// WIDTH (CLHEP::mm): diameter // REPEAT: Repeat this many times in phi with equal distance between them. // BOX // Ignored: RIN2, ROUT2 @@ -93,9 +93,8 @@ #include <iomanip> namespace InDetDD { - -ServiceVolume::ServiceVolume() - : m_rmin(0), + ServiceVolume::ServiceVolume() + : m_rmin(0), m_rmax(0), m_rmin2(0), m_rmax2(0), @@ -120,272 +119,252 @@ ServiceVolume::ServiceVolume() m_envNum(0), m_envParentNum(0), m_zShift(0.) -{} + {} -void -ServiceVolume::reduceSize(double safety) -{ - // Don't do anything if its a very thin volume. - if (length() > 4.*safety) { - if (m_zmax < m_zmin) std::swap(m_zmin,m_zmax); - m_safety = safety; + void + ServiceVolume::reduceSize(double safety) { + // Don't do anything if its a very thin volume. + if (length() > 4. * safety) { + if (m_zmax < m_zmin) std::swap(m_zmin, m_zmax); + m_safety = safety; + } + m_geoShape = 0; } - m_geoShape = 0; -} - -void -ServiceVolume::setLabel(const std::string & name, int volId) -{ - std::ostringstream o; - o.fill('0'); - o << name << std::setw(2) << volId; - m_label = o.str(); -} - -std::string -ServiceVolume::fullLabel() const -{ - if (m_volName.empty()) return m_label; - return m_label+"_"+m_volName; -} -void -ServiceVolume::print() const -{ - std::cout << m_rmin << " " - << m_rmax << " " - << m_zmin << " " - << m_zmax << " " - << m_region << " " - << fullLabel() - << std::endl; -} - -const GeoShape * -ServiceVolume::getShape() const { - - // If prebuilt then return - if (m_geoShape.get()) return m_geoShape.get(); - - // - // Dimensions - // - //double rmin = rmin(); - //double rmax = rmax(); - //double rmin2 = rmin2(); - //double rmax2 = rmax2(); - //double phiLoc = phiLoc(); - //double phiWidth = phiWidth(); - //int sides = sides(); - //const std::string & shapeType = shapeType(); + void + ServiceVolume::setLabel(const std::string& name, int volId) { + std::ostringstream o; + o.fill('0'); + o << name << std::setw(2) << volId; + m_label = o.str(); + } - double halflength = 0.5 * length(); + std::string + ServiceVolume::fullLabel() const { + if (m_volName.empty()) return m_label; - //std::cout << "Building service volume " << logName << ": " - // << rmin << ", " - // << rmax << ", " - // << halflength << ", " - // << materialName << std::endl; + return m_label + "_" + m_volName; + } - const GeoShape * serviceShape = 0; - double volume = 0; + void + ServiceVolume::print() const { + std::cout << m_rmin << " " + << m_rmax << " " + << m_zmin << " " + << m_zmax << " " + << m_region << " " + << fullLabel() + << std::endl; + } - // Check if service needs to be shifted - // if(fabs(m_zShift)>0.001) - // std::cout<<"SHIFTED SERVICE : "<<m_volName<<" "<<m_shapeType<<std::endl; - - if (m_shapeType.empty() || m_shapeType == "TUBE") { - serviceShape = new GeoTube(m_rmin,m_rmax,halflength); - } else if (m_shapeType == "TUBS") { - serviceShape = new GeoTubs(m_rmin,m_rmax,halflength,m_phiLoc,m_phiWidth); - } else if (m_shapeType == "CONS" || m_shapeType == "CONE") { - double phiWidthTmp = m_phiWidth; - if (m_shapeType == "CONE" || phiWidthTmp == 0) { - phiWidthTmp = 2 * M_PI; - } - serviceShape = new GeoCons(m_rmin,m_rmin2,m_rmax,m_rmax2,halflength,m_phiLoc,phiWidthTmp); - } else if (m_shapeType == "PGON") { - GeoPgon * shapeTmp = new GeoPgon(m_phiLoc,2*M_PI,m_sides); - shapeTmp->addPlane(-halflength,m_rmin,m_rmax); - shapeTmp->addPlane(halflength,m_rmin2,m_rmax2); - serviceShape = shapeTmp; - } else if (m_shapeType == "PGON2") { - // Radius defined at the side, not the corner - double alpha = M_PI/m_sides; - double cosalpha = cos(alpha); - double rminB = m_rmin/cosalpha; - double rmaxB = m_rmax/cosalpha; - double rmin2B = m_rmin2/cosalpha; - double rmax2B = m_rmax2/cosalpha; - GeoPgon * shapeTmp = new GeoPgon(m_phiLoc-alpha,2*M_PI,m_sides); - shapeTmp->addPlane(-halflength,rminB,rmaxB); - shapeTmp->addPlane(halflength,rmin2B,rmax2B); - serviceShape = shapeTmp; - } else if (m_shapeType == "PGON3" || m_shapeType == "PGON4") { - // Outer edge - GeoPgon * shapeTmp1 = 0; - if (m_shapeType == "PGON3") { - shapeTmp1 = new GeoPgon(m_phiLoc,2*M_PI,m_sides); - shapeTmp1->addPlane(-halflength,0,m_rmax); - shapeTmp1->addPlane(halflength,0,m_rmax2); - } else { //PGON4 - double alpha = M_PI/m_sides; - double cosalpha = cos(alpha); - double rmaxB = m_rmax/cosalpha; - double rmax2B = m_rmax2/cosalpha; - shapeTmp1 = new GeoPgon(m_phiLoc-alpha,2*M_PI,m_sides); - shapeTmp1->addPlane(-halflength,0,rmaxB); - shapeTmp1->addPlane(halflength,0,rmax2B); + const GeoShape* + ServiceVolume::getShape() const { + // If prebuilt then return + if (m_geoShape.get()) return m_geoShape.get(); + + // + // Dimensions + // + //double rmin = rmin(); + //double rmax = rmax(); + //double rmin2 = rmin2(); + //double rmax2 = rmax2(); + //double phiLoc = phiLoc(); + //double phiWidth = phiWidth(); + //int sides = sides(); + //const std::string & shapeType = shapeType(); + + double halflength = 0.5 * length(); + + //std::cout << "Building service volume " << logName << ": " + // << rmin << ", " + // << rmax << ", " + // << halflength << ", " + // << materialName << std::endl; + + const GeoShape* serviceShape = 0; + double volume = 0; + + // Check if service needs to be shifted + // if(fabs(m_zShift)>0.001) + // std::cout<<"SHIFTED SERVICE : "<<m_volName<<" "<<m_shapeType<<std::endl; + + if (m_shapeType.empty() || m_shapeType == "TUBE") { + serviceShape = new GeoTube(m_rmin, m_rmax, halflength); + } else if (m_shapeType == "TUBS") { + serviceShape = new GeoTubs(m_rmin, m_rmax, halflength, m_phiLoc, m_phiWidth); + } else if (m_shapeType == "CONS" || m_shapeType == "CONE") { + double phiWidthTmp = m_phiWidth; + if (m_shapeType == "CONE" || phiWidthTmp == 0) { + phiWidthTmp = 2 * M_PI; + } + serviceShape = new GeoCons(m_rmin, m_rmin2, m_rmax, m_rmax2, halflength, m_phiLoc, phiWidthTmp); + } else if (m_shapeType == "PGON") { + GeoPgon* shapeTmp = new GeoPgon(m_phiLoc, 2 * M_PI, m_sides); + shapeTmp->addPlane(-halflength, m_rmin, m_rmax); + shapeTmp->addPlane(halflength, m_rmin2, m_rmax2); + serviceShape = shapeTmp; + } else if (m_shapeType == "PGON2") { + // Radius defined at the side, not the corner + double alpha = M_PI / m_sides; + double cosalpha = cos(alpha); + double rminB = m_rmin / cosalpha; + double rmaxB = m_rmax / cosalpha; + double rmin2B = m_rmin2 / cosalpha; + double rmax2B = m_rmax2 / cosalpha; + GeoPgon* shapeTmp = new GeoPgon(m_phiLoc - alpha, 2 * M_PI, m_sides); + shapeTmp->addPlane(-halflength, rminB, rmaxB); + shapeTmp->addPlane(halflength, rmin2B, rmax2B); + serviceShape = shapeTmp; + } else if (m_shapeType == "PGON3" || m_shapeType == "PGON4") { + // Outer edge + GeoPgon* shapeTmp1 = 0; + if (m_shapeType == "PGON3") { + shapeTmp1 = new GeoPgon(m_phiLoc, 2 * M_PI, m_sides); + shapeTmp1->addPlane(-halflength, 0, m_rmax); + shapeTmp1->addPlane(halflength, 0, m_rmax2); + } else { //PGON4 + double alpha = M_PI / m_sides; + double cosalpha = cos(alpha); + double rmaxB = m_rmax / cosalpha; + double rmax2B = m_rmax2 / cosalpha; + shapeTmp1 = new GeoPgon(m_phiLoc - alpha, 2 * M_PI, m_sides); + shapeTmp1->addPlane(-halflength, 0, rmaxB); + shapeTmp1->addPlane(halflength, 0, rmax2B); + } + // Don't trust boolean volume calculation. + volume = shapeTmp1->volume(); + // Inner edge + GeoShape* shapeTmp2 = 0; + if (m_rmin == m_rmin2) { + shapeTmp2 = new GeoTube(0, m_rmin, halflength + 0.1 * CLHEP::mm); + volume -= 2 * M_PI * m_rmin * m_rmin * halflength; + } else { + shapeTmp2 = new GeoCons(0, 0, m_rmin, m_rmin2, halflength + 0.1 * CLHEP::mm, 0, 2 * M_PI); + volume -= 2 * M_PI * pow(0.5 * (m_rmin + m_rmin2), 2) * halflength; + } + serviceShape = &(shapeTmp1->subtract(*shapeTmp2)); } - // Don't trust boolean volume calculation. - volume = shapeTmp1->volume(); - // Inner edge - GeoShape * shapeTmp2 = 0; - if (m_rmin == m_rmin2) { - shapeTmp2 = new GeoTube(0,m_rmin,halflength+0.1*CLHEP::mm); - volume -= 2*M_PI*m_rmin*m_rmin*halflength; - } else { - shapeTmp2 = new GeoCons(0,0,m_rmin,m_rmin2,halflength+0.1*CLHEP::mm,0,2*M_PI); - volume -= 2*M_PI*pow(0.5*(m_rmin+m_rmin2),2)*halflength; - } - serviceShape = &(shapeTmp1->subtract(*shapeTmp2)); - } // else if (m_shapeType == "PGON31"){ // // Outer edge // GeoTube *shapeTmp1 = new GeoTube(0,m_rmax,halflength); - // halflength+=0.1*CLHEP::mm; // double alpha = M_PI/m_sides; -// double cosalpha = cos(alpha); +// double cosalpha = cos(alpha); // double rmaxB = m_rmin/cosalpha; // double rmax2B = m_rmin2/cosalpha; // GeoPgon* shapeTmp2 = new GeoPgon(m_phiLoc-alpha,2*M_PI,m_sides); // shapeTmp2->addPlane(-halflength,0.,rmaxB); // shapeTmp2->addPlane(halflength,0.,rmax2B); - -// // Don't trust boolean volume calculation. +// // Don't trust boolean volume calculation. // volume = shapeTmp1->volume() - shapeTmp2->volume(); // serviceShape = &(shapeTmp1->subtract(*shapeTmp2)); -// } - else if (m_shapeType == "ROD") { - serviceShape = new GeoTube(0,0.5*m_phiWidth,halflength); - } - else if (m_shapeType == "ROD2") { - // std::cout<<"ROD2 : "<<m_rmin<<" "<<m_rmin2<<" "<<0.5*m_phiWidth<<" "<<halflength<<std::endl; - serviceShape = new GeoTube(m_rmin2-m_rmin,0.5*m_phiWidth,halflength); - } - else if (m_shapeType == "BOX") { - serviceShape = new GeoBox(0.5*(m_rmax-m_rmin),0.5*m_phiWidth,halflength); - } - else if (m_shapeType == "TRAP") { - double thickness = 0.5*(m_rmax-m_rmin); - double averad = 0.5*(m_rmin+m_rmax); - double w1 = 0.5*m_phiWidth*m_rmin/averad; - double w2 = 0.5*m_phiWidth*m_rmax/averad; - serviceShape = new GeoTrap(halflength, 0, 0, thickness, w1, w2, 0, thickness, w1, w2, 0); - } else { - // msg(MSG::ERROR) << "Unrecognized shape for services" << m_shapeType << endmsg; - std::cout << "ServiceVolume: ERROR: Unrecognized shape for services" << m_shapeType << std::endl; - } - - if (!volume&&serviceShape!=0) volume = serviceShape->volume(); - - m_volume = volume; - m_geoShape.set(serviceShape); - return serviceShape; -} +// } + else if (m_shapeType == "ROD") { + serviceShape = new GeoTube(0, 0.5 * m_phiWidth, halflength); + } else if (m_shapeType == "ROD2") { + // std::cout<<"ROD2 : "<<m_rmin<<" "<<m_rmin2<<" "<<0.5*m_phiWidth<<" "<<halflength<<std::endl; + serviceShape = new GeoTube(m_rmin2 - m_rmin, 0.5 * m_phiWidth, halflength); + } else if (m_shapeType == "BOX") { + serviceShape = new GeoBox(0.5 * (m_rmax - m_rmin), 0.5 * m_phiWidth, halflength); + } else if (m_shapeType == "TRAP") { + double thickness = 0.5 * (m_rmax - m_rmin); + double averad = 0.5 * (m_rmin + m_rmax); + double w1 = 0.5 * m_phiWidth * m_rmin / averad; + double w2 = 0.5 * m_phiWidth * m_rmax / averad; + serviceShape = new GeoTrap(halflength, 0, 0, thickness, w1, w2, 0, thickness, w1, w2, 0); + } else { + // msg(MSG::ERROR) << "Unrecognized shape for services" << m_shapeType << endmsg; + std::cout << "ServiceVolume: ERROR: Unrecognized shape for services" << m_shapeType << std::endl; + } -double -ServiceVolume::volume() const -{ - // Make sure shape is already built. - getShape(); - return m_volume; -} + if (!volume && serviceShape != 0) volume = serviceShape->volume(); -void -ServiceVolume::setGeoShape(const GeoShape * geoShape, double volume) -{ - m_geoShape.reset(); - if (geoShape) { m_volume = volume; - // We allow a volume to specified as the volume calculation for some shapes (ie boolean volumes) are unreliable. - // If volume is not supplied, get it from the shape itself. - if (!m_volume) m_volume = geoShape->volume(); - m_geoShape.set(geoShape); - m_lockGeoShape = true; // This disables resetGeoShape(). - setShapeType("CUSTOM"); - } else { - // If pass null pointer we unlock the shape. - m_lockGeoShape = false; + m_geoShape.set(serviceShape); + return serviceShape; } -} + double + ServiceVolume::volume() const { + // Make sure shape is already built. + getShape(); + return m_volume; + } -double -ServiceVolume::origVolume() const -{ - if (m_origVolume) return m_origVolume; - return volume(); -} + void + ServiceVolume::setGeoShape(const GeoShape* geoShape, double volume) { + m_geoShape.reset(); + if (geoShape) { + m_volume = volume; + // We allow a volume to specified as the volume calculation for some shapes (ie boolean volumes) are unreliable. + // If volume is not supplied, get it from the shape itself. + if (!m_volume) m_volume = geoShape->volume(); + m_geoShape.set(geoShape); + m_lockGeoShape = true; // This disables resetGeoShape(). + setShapeType("CUSTOM"); + } else { + // If pass null pointer we unlock the shape. + m_lockGeoShape = false; + } + } -void -ServiceVolume::setSplittable() -{ - m_splittableR = m_splittableZ = true; - if (m_shapeType == "CUSTOM") { - m_splittableR = m_splittableZ = false; - } else if (!(m_shapeType == "" || m_shapeType == "TUBE" || m_shapeType == "TUBS")) { - m_splittableR = false; - } -} + double + ServiceVolume::origVolume() const { + if (m_origVolume) return m_origVolume; -GeoShapeHolder::GeoShapeHolder() - : m_geoShape(0) -{} + return volume(); + } -GeoShapeHolder::GeoShapeHolder(const GeoShape * geoShape) - : m_geoShape(geoShape) -{ - if (m_geoShape) m_geoShape->ref(); -} + void + ServiceVolume::setSplittable() { + m_splittableR = m_splittableZ = true; + if (m_shapeType == "CUSTOM") { + m_splittableR = m_splittableZ = false; + } else if (!(m_shapeType == "" || m_shapeType == "TUBE" || m_shapeType == "TUBS")) { + m_splittableR = false; + } + } -GeoShapeHolder::~GeoShapeHolder() -{ - reset(); -} + GeoShapeHolder::GeoShapeHolder() + : m_geoShape(0) + {} - -GeoShapeHolder::GeoShapeHolder(const GeoShapeHolder & rhs) -{ - m_geoShape = rhs.m_geoShape; - if (m_geoShape) m_geoShape->ref(); -} + GeoShapeHolder::GeoShapeHolder(const GeoShape* geoShape) + : m_geoShape(geoShape) { + if (m_geoShape) m_geoShape->ref(); + } -GeoShapeHolder & GeoShapeHolder::operator=(const GeoShapeHolder & rhs) -{ - if (&rhs != this) { + GeoShapeHolder::~GeoShapeHolder() { reset(); + } + + GeoShapeHolder::GeoShapeHolder(const GeoShapeHolder& rhs) { m_geoShape = rhs.m_geoShape; if (m_geoShape) m_geoShape->ref(); } - return *this; -} - -void GeoShapeHolder::set(const GeoShape * geoShape) -{ - reset(); - m_geoShape = geoShape; - if (m_geoShape) m_geoShape->ref(); -} -void GeoShapeHolder::reset() -{ - if (m_geoShape) m_geoShape->unref(); - m_geoShape = 0; -} + GeoShapeHolder& + GeoShapeHolder::operator = (const GeoShapeHolder& rhs) { + if (&rhs != this) { + reset(); + m_geoShape = rhs.m_geoShape; + if (m_geoShape) m_geoShape->ref(); + } + return *this; + } + void + GeoShapeHolder::set(const GeoShape* geoShape) { + reset(); + m_geoShape = geoShape; + if (m_geoShape) m_geoShape->ref(); + } + void + GeoShapeHolder::reset() { + if (m_geoShape) m_geoShape->unref(); + m_geoShape = 0; + } } // end namespace - diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ServiceVolumeMaker.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ServiceVolumeMaker.cxx old mode 100644 new mode 100755 index 95cbc03c408a95530e75b7ab4e5c5130c044862e..c980d5cc37dc21ae2661a1da58a2401bd8dd6352 --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ServiceVolumeMaker.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/ServiceVolumeMaker.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/ServiceVolumeMaker.h" #include "InDetGeoModelUtils/ServiceVolume.h" @@ -13,402 +13,391 @@ #include "CLHEP/Units/SystemOfUnits.h" namespace InDetDD { + ServiceVolumeSchema::ServiceVolumeSchema() { + setSimpleSchema(); + } + + void + ServiceVolumeSchema::setPixelSchema() { + m_rmin = "RIN"; + m_rmax = "ROUT"; + m_rmin2 = "RIN2"; + m_rmax2 = "ROUT2"; + m_zmin = "ZIN"; + m_zmax = "ZOUT"; + m_zsymm = "ZSYMM"; + m_materialName = "MATERIALNAME"; + m_repeat = "REPEAT"; + m_phiStart = "PHI"; + m_phiDelta = "WIDTH"; + m_width = "WIDTH"; + m_shapeType = "SHAPE"; + m_volName = "VOLNAME"; + m_radialDiv = ""; + m_phiStep = ""; + m_volId = "FRAMENUM"; + m_shiftFlag = "SHIFT"; + } + + void + ServiceVolumeSchema::setDefaultSchema() { + m_rmin = "RMIN"; + m_rmax = "RMAX"; + m_rmin2 = "RMIN2"; + m_rmax2 = "RMAX2"; + m_zmin = "ZMIN"; + m_zmax = "ZMAX"; + m_zsymm = "ZSYMM"; + m_materialName = "MATERIAL"; + m_repeat = "NREPEAT"; + m_phiStart = "PHISTART"; + m_phiDelta = "PHIDELTA"; + m_width = ""; + m_shapeType = ""; + m_volName = "NAME"; + m_radialDiv = "RADIAL"; + m_phiStep = "PHISTEP"; + m_volId = ""; + m_shiftFlag = "SHIFT"; + } + + void + ServiceVolumeSchema::setSimpleSchema() { + m_rmin = "RMIN"; + m_rmax = "RMAX"; + m_rmin2 = ""; + m_rmax2 = ""; + m_zmin = "ZMIN"; + m_zmax = "ZMAX"; + m_zsymm = "ZSYMM"; + m_materialName = "MATERIAL"; + m_repeat = ""; + m_phiStart = ""; + m_phiDelta = ""; + m_width = ""; + m_shapeType = ""; + m_volName = "NAME"; + m_radialDiv = ""; + m_phiStep = ""; + m_volId = ""; + m_shiftFlag = "SHIFT"; + } -ServiceVolumeSchema::ServiceVolumeSchema() -{ - setSimpleSchema(); -} - -void ServiceVolumeSchema::setPixelSchema() -{ - m_rmin = "RIN"; - m_rmax = "ROUT"; - m_rmin2 = "RIN2"; - m_rmax2 = "ROUT2"; - m_zmin = "ZIN"; - m_zmax = "ZOUT"; - m_zsymm = "ZSYMM"; - m_materialName = "MATERIALNAME"; - m_repeat = "REPEAT"; - m_phiStart = "PHI"; - m_phiDelta = "WIDTH"; - m_width = "WIDTH"; - m_shapeType = "SHAPE"; - m_volName = "VOLNAME"; - m_radialDiv = ""; - m_phiStep = ""; - m_volId = "FRAMENUM"; - m_shiftFlag="SHIFT"; -} - -void ServiceVolumeSchema::setDefaultSchema() -{ - m_rmin = "RMIN"; - m_rmax = "RMAX"; - m_rmin2 = "RMIN2"; - m_rmax2 = "RMAX2"; - m_zmin = "ZMIN"; - m_zmax = "ZMAX"; - m_zsymm = "ZSYMM"; - m_materialName = "MATERIAL"; - m_repeat = "NREPEAT"; - m_phiStart = "PHISTART"; - m_phiDelta = "PHIDELTA"; - m_width = ""; - m_shapeType = ""; - m_volName = "NAME"; - m_radialDiv = "RADIAL"; - m_phiStep = "PHISTEP"; - m_volId = ""; - m_shiftFlag="SHIFT"; -} - -void ServiceVolumeSchema::setSimpleSchema() -{ - m_rmin = "RMIN"; - m_rmax = "RMAX"; - m_rmin2 = ""; - m_rmax2 = ""; - m_zmin = "ZMIN"; - m_zmax = "ZMAX"; - m_zsymm = "ZSYMM"; - m_materialName = "MATERIAL"; - m_repeat = ""; - m_phiStart = ""; - m_phiDelta = ""; - m_width = ""; - m_shapeType = ""; - m_volName = "NAME"; - m_radialDiv = ""; - m_phiStep = ""; - m_volId = ""; - m_shiftFlag="SHIFT"; - -} - -ServiceVolumeMakerMgr::ServiceVolumeMakerMgr(IRDBRecordset_ptr table, const ServiceVolumeSchema & schema, - const InDetDD::AthenaComps * athenaComps) - : m_table(table), + ServiceVolumeMakerMgr::ServiceVolumeMakerMgr(IRDBRecordset_ptr table, const ServiceVolumeSchema& schema, + const InDetDD::AthenaComps* athenaComps) + : m_table(table), m_schema(schema), m_athenaComps(athenaComps) -{} - -const IGeometryDBSvc * -ServiceVolumeMakerMgr::db() const { - return m_athenaComps->geomDB(); -} - -double ServiceVolumeMakerMgr::rmin(int index) const -{ - return db()->getDouble(m_table, m_schema.rmin(), index) * CLHEP::mm; -} - - -double ServiceVolumeMakerMgr::rmax(int index) const -{ - return db()->getDouble(m_table, m_schema.rmax(), index) * CLHEP::mm; -} - - -double ServiceVolumeMakerMgr::rmin2(int index) const -{ - return db()->getDouble(m_table, m_schema.rmin2(), index) * CLHEP::mm; -} - -double ServiceVolumeMakerMgr::rmax2(int index) const -{ - return db()->getDouble(m_table, m_schema.rmax2(), index) * CLHEP::mm; -} - -double ServiceVolumeMakerMgr::zmin(int index) const -{ - return db()->getDouble(m_table, m_schema.zmin(), index) * CLHEP::mm; -} - -double ServiceVolumeMakerMgr::zmax(int index) const -{ - return db()->getDouble(m_table, m_schema.zmax(), index) * CLHEP::mm; -} - -double ServiceVolumeMakerMgr::phiDelta(int index) const -{ - return db()->getDouble(m_table, m_schema.phiDelta(), index) * CLHEP::deg; -} - -double ServiceVolumeMakerMgr::width(int index) const -{ - if (m_schema.has_width()) { - return db()->getDouble(m_table, m_schema.width(), index) * CLHEP::mm; + {} + + const IGeometryDBSvc* + ServiceVolumeMakerMgr::db() const { + return m_athenaComps->geomDB(); + } + + double + ServiceVolumeMakerMgr::rmin(int index) const { + return db()->getDouble(m_table, m_schema.rmin(), index) * CLHEP::mm; + } + + double + ServiceVolumeMakerMgr::rmax(int index) const { + return db()->getDouble(m_table, m_schema.rmax(), index) * CLHEP::mm; + } + + double + ServiceVolumeMakerMgr::rmin2(int index) const { + return db()->getDouble(m_table, m_schema.rmin2(), index) * CLHEP::mm; + } + + double + ServiceVolumeMakerMgr::rmax2(int index) const { + return db()->getDouble(m_table, m_schema.rmax2(), index) * CLHEP::mm; + } + + double + ServiceVolumeMakerMgr::zmin(int index) const { + return db()->getDouble(m_table, m_schema.zmin(), index) * CLHEP::mm; + } + + double + ServiceVolumeMakerMgr::zmax(int index) const { + return db()->getDouble(m_table, m_schema.zmax(), index) * CLHEP::mm; + } + + double + ServiceVolumeMakerMgr::phiDelta(int index) const { + return db()->getDouble(m_table, m_schema.phiDelta(), index) * CLHEP::deg; + } + + double + ServiceVolumeMakerMgr::width(int index) const { + if (m_schema.has_width()) { + return db()->getDouble(m_table, m_schema.width(), index) * CLHEP::mm; + } + return 0; + } + + double + ServiceVolumeMakerMgr::phiStart(int index) const { + return db()->getDouble(m_table, m_schema.phiStart(), index) * CLHEP::deg; } - return 0; -} - -double ServiceVolumeMakerMgr::phiStart(int index) const -{ - return db()->getDouble(m_table, m_schema.phiStart(), index) * CLHEP::deg; -} - -double ServiceVolumeMakerMgr::phiStep(int index) const -{ - if (m_schema.has_phiStep()) { - return db()->getDouble(m_table, m_schema.phiStep(), index) * CLHEP::deg; - } - return 0; -} - -bool ServiceVolumeMakerMgr::zsymm(int index) const -{ - return db()->getInt(m_table, m_schema.zsymm(), index); -} - - -int ServiceVolumeMakerMgr::repeat(int index) const -{ - return db()->getInt(m_table, m_schema.repeat(), index); -} - -int ServiceVolumeMakerMgr::radialDiv(int index) const -{ - if (m_schema.has_radial()) { - return db()->getInt(m_table, m_schema.radialDiv(), index); - } else { + + double + ServiceVolumeMakerMgr::phiStep(int index) const { + if (m_schema.has_phiStep()) { + return db()->getDouble(m_table, m_schema.phiStep(), index) * CLHEP::deg; + } return 0; } -} -std::string ServiceVolumeMakerMgr::shapeType(int index) const -{ - if (m_schema.has_shapeType()) { - if (db()->testField(m_table, m_schema.shapeType(), index)) { - return db()->getString(m_table, m_schema.shapeType(), index); + bool + ServiceVolumeMakerMgr::zsymm(int index) const { + return db()->getInt(m_table, m_schema.zsymm(), index); + } + + int + ServiceVolumeMakerMgr::repeat(int index) const { + return db()->getInt(m_table, m_schema.repeat(), index); + } + + int + ServiceVolumeMakerMgr::radialDiv(int index) const { + if (m_schema.has_radial()) { + return db()->getInt(m_table, m_schema.radialDiv(), index); } else { - return "TUBE"; + return 0; + } + } + + std::string + ServiceVolumeMakerMgr::shapeType(int index) const { + if (m_schema.has_shapeType()) { + if (db()->testField(m_table, m_schema.shapeType(), index)) { + return db()->getString(m_table, m_schema.shapeType(), index); + } else { + return "TUBE"; + } + } + return "UNKNOWN"; + } + + std::string + ServiceVolumeMakerMgr::volName(int index) const { + if (db()->testField(m_table, m_schema.volName(), index)) { + return db()->getString(m_table, m_schema.volName(), index); } + return ""; } - return "UNKNOWN"; -} - -std::string ServiceVolumeMakerMgr::volName(int index) const -{ - if (db()->testField(m_table, m_schema.volName(), index)) { - return db()->getString(m_table, m_schema.volName(), index); - } - return ""; -} - -std::string ServiceVolumeMakerMgr::materialName(int index) const -{ - return db()->getString(m_table, m_schema.materialName(), index); -} - -unsigned int ServiceVolumeMakerMgr::numElements() const -{ - return db()->getTableSize(m_table); -} - -int ServiceVolumeMakerMgr::volId(int index) const -{ - if (m_schema.has_volId()) { - return db()->getInt(m_table, m_schema.volId(), index); + + std::string + ServiceVolumeMakerMgr::materialName(int index) const { + return db()->getString(m_table, m_schema.materialName(), index); + } + + unsigned int + ServiceVolumeMakerMgr::numElements() const { + return db()->getTableSize(m_table); + } + + int + ServiceVolumeMakerMgr::volId(int index) const { + if (m_schema.has_volId()) { + return db()->getInt(m_table, m_schema.volId(), index); + } + return 0; } - return 0; -} - -int ServiceVolumeMakerMgr::shiftFlag(int index) const -{ - if (m_schema.has_shiftFlag()) { - if (db()->testField(m_table, m_schema.shiftFlag(), index)) - return db()->getInt(m_table, m_schema.shiftFlag(), index); + + int + ServiceVolumeMakerMgr::shiftFlag(int index) const { + if (m_schema.has_shiftFlag()) { + if (db()->testField(m_table, m_schema.shiftFlag(), index)) return db()->getInt(m_table, + m_schema.shiftFlag(), index); + } + return 0; } - return 0; -} -std::vector<double> ServiceVolumeMakerMgr::readLayerShift() const -{ - std::vector<double> layerShift; + std::vector<double> + ServiceVolumeMakerMgr::readLayerShift() const { + std::vector<double> layerShift; - IRDBAccessSvc *rdbSvc = m_athenaComps->rdbAccessSvc(); - IGeoDbTagSvc *geoDbTag = m_athenaComps->geoDbTagSvc(); + IRDBAccessSvc* rdbSvc = m_athenaComps->rdbAccessSvc(); + IGeoDbTagSvc* geoDbTag = m_athenaComps->geoDbTagSvc(); - DecodeVersionKey versionKey(geoDbTag,"Pixel"); - std::string detectorKey = versionKey.tag(); - std::string detectorNode = versionKey.node(); + DecodeVersionKey versionKey(geoDbTag, "Pixel"); + std::string detectorKey = versionKey.tag(); + std::string detectorNode = versionKey.node(); - IRDBRecordset_ptr PixelBarrelGeneral= rdbSvc->getRecordsetPtr("PixelBarrelGeneral",detectorKey, detectorNode); - IRDBRecordset_ptr PixelLayer= rdbSvc->getRecordsetPtr("PixelLayer",detectorKey, detectorNode); + IRDBRecordset_ptr PixelBarrelGeneral = rdbSvc->getRecordsetPtr("PixelBarrelGeneral", detectorKey, detectorNode); + IRDBRecordset_ptr PixelLayer = rdbSvc->getRecordsetPtr("PixelLayer", detectorKey, detectorNode); - int numLayers = db()->getInt(PixelBarrelGeneral,"NLAYER"); - for(int iLayer=0; iLayer<numLayers; iLayer++) - { + int numLayers = db()->getInt(PixelBarrelGeneral, "NLAYER"); + for (int iLayer = 0; iLayer < numLayers; iLayer++) { double shift = 0; - if (db()->testField(PixelLayer,"GBLSHIFT",iLayer)) shift = db()->getDouble(PixelLayer,"GBLSHIFT",iLayer); + if (db()->testField(PixelLayer, "GBLSHIFT", iLayer)) shift = db()->getDouble(PixelLayer, "GBLSHIFT", iLayer); layerShift.push_back(shift); } - - return layerShift; -} - -ServiceVolumeMaker::ServiceVolumeMaker(const std::string & label, - IRDBRecordset_ptr table, const ServiceVolumeSchema & schema, - const InDetDD::AthenaComps * athenaComps) - : m_label(label) -{ - m_mgr = new ServiceVolumeMakerMgr(table, schema, athenaComps); - m_layerShift = m_mgr->readLayerShift(); - // std::cout<<"LAYER SHIFT "<<m_layerShift[0]<<" "<<m_layerShift[1]<<" "<<m_layerShift[2]<<" "<<m_layerShift[3]<<std::endl; -} - -ServiceVolumeMaker::~ServiceVolumeMaker() -{ - for (unsigned int i = 0; i < m_services.size(); i++) { - delete m_services[i]; - } - delete m_mgr; -} - -const std::vector<const ServiceVolume *> & -ServiceVolumeMaker::makeAll() -{ - for (unsigned int ii = 0; ii < numElements(); ++ii) { - m_services.push_back(make(ii)); + + return layerShift; } - return m_services; -} - -unsigned int -ServiceVolumeMaker::numElements() const { - return m_mgr->numElements(); -} - -ServiceVolume * -ServiceVolumeMaker::make(int ii) -{ - // - // Retrieve/calculate the parameters for the volume. - // - ServiceVolume * param = new ServiceVolume ; - param->setMaterial(m_mgr->materialName(ii)); - param->setRmin(m_mgr->rmin(ii)); - param->setRmax(m_mgr->rmax(ii)); - param->setZmin(m_mgr->zmin(ii)); - param->setZmax(m_mgr->zmax(ii)); - param->setZsymm(m_mgr->zsymm(ii)); - param->setVolName(m_mgr->volName(ii)); - - double zShift=0.; // the famous IBL Z shift - if(m_mgr->shiftFlag(ii)>0) zShift=m_layerShift[m_mgr->shiftFlag(ii)-100]; - param->setZShift(zShift); - - int volId = m_mgr->volId(ii); - if (volId == 0) volId = ii+1; - - bool needsRotation = false; - - // For TUBE there is no need to read the rest - std::string shapeType = m_mgr->shapeType(ii); - if (!m_mgr->schema().simple() && !shapeType.empty() && shapeType != "TUBE") { - - double rmin2 = m_mgr->rmin2(ii); - double rmax2 = m_mgr->rmax2(ii); - - if (rmin2 <= 0) rmin2 = param->rmin(); - if (rmax2 <= 0) rmax2 = param->rmax(); - - int radialDiv = m_mgr->radialDiv(ii); - - double phiDelta = m_mgr->phiDelta(ii); - - bool fullPhiSector = false; - if (phiDelta == 0 || phiDelta >=359.9*CLHEP::degree) { - phiDelta = 360*CLHEP::degree; - fullPhiSector = true; - } - //else { - //phiDelta -= 2*phiepsilon; - //phiStart += phiepsilon; - // } - - if (shapeType == "UNKNOWN") { - if (radialDiv > 0) { - shapeType = "RADIAL"; - } else if (param->rmin() == rmin2 && param->rmax() == rmax2 ) { - if (fullPhiSector) { - shapeType = "TUBE"; - } else { - shapeType = "TUBS"; + + ServiceVolumeMaker::ServiceVolumeMaker(const std::string& label, + IRDBRecordset_ptr table, const ServiceVolumeSchema& schema, + const InDetDD::AthenaComps* athenaComps) + : m_label(label) { + m_mgr = new ServiceVolumeMakerMgr(table, schema, athenaComps); + m_layerShift = m_mgr->readLayerShift(); + // std::cout<<"LAYER SHIFT "<<m_layerShift[0]<<" "<<m_layerShift[1]<<" "<<m_layerShift[2]<<" + // "<<m_layerShift[3]<<std::endl; } - } else { - shapeType = "CONS"; - } + + ServiceVolumeMaker::~ServiceVolumeMaker() { + for (unsigned int i = 0; i < m_services.size(); i++) { + delete m_services[i]; } - - - int repeat = m_mgr->repeat(ii); - if (repeat == 0) repeat = 1; - - double phiStart = m_mgr->phiStart(ii); - double phiWidth = phiDelta; - - if (shapeType == "CONS" || shapeType == "TUBS") { - const double phiepsilon = 0.001*CLHEP::degree; - phiWidth -= 2*phiepsilon; - phiStart += phiepsilon; + delete m_mgr; + } + + const std::vector<const ServiceVolume*>& + ServiceVolumeMaker::makeAll() { + for (unsigned int ii = 0; ii < numElements(); ++ii) { + m_services.push_back(make(ii)); } - - // Can be in degree or CLHEP::mm. Usually it is CLHEP::deg expect for BOX, TRAP and ROD shape - // Geometry manager makes no assumptions about units. So we must interpret here. - if (shapeType == "BOX" || shapeType == "ROD" || shapeType=="ROD2" || shapeType == "TRAP") { - phiWidth = m_mgr->width(ii); // in mm - } - - if (shapeType == "PGON" || shapeType == "PGON2" || - shapeType == "CONE" || shapeType == "CONS" || - shapeType == "PGON3" || shapeType == "PGON4") { - if ((rmin2 != param->rmin()) || (rmax2 != param->rmax())) { - needsRotation = true; + return m_services; + } + + unsigned int + ServiceVolumeMaker::numElements() const { + return m_mgr->numElements(); + } + + ServiceVolume* + ServiceVolumeMaker::make(int ii) { + // + // Retrieve/calculate the parameters for the volume. + // + ServiceVolume* param = new ServiceVolume; + + param->setMaterial(m_mgr->materialName(ii)); + param->setRmin(m_mgr->rmin(ii)); + param->setRmax(m_mgr->rmax(ii)); + param->setZmin(m_mgr->zmin(ii)); + param->setZmax(m_mgr->zmax(ii)); + param->setZsymm(m_mgr->zsymm(ii)); + param->setVolName(m_mgr->volName(ii)); + + double zShift = 0.; // the famous IBL Z shift + if (m_mgr->shiftFlag(ii) > 0) zShift = m_layerShift[m_mgr->shiftFlag(ii) - 100]; + param->setZShift(zShift); + + int volId = m_mgr->volId(ii); + if (volId == 0) volId = ii + 1; + + bool needsRotation = false; + + // For TUBE there is no need to read the rest + std::string shapeType = m_mgr->shapeType(ii); + if (!m_mgr->schema().simple() && !shapeType.empty() && shapeType != "TUBE") { + double rmin2 = m_mgr->rmin2(ii); + double rmax2 = m_mgr->rmax2(ii); + + if (rmin2 <= 0) rmin2 = param->rmin(); + if (rmax2 <= 0) rmax2 = param->rmax(); + + int radialDiv = m_mgr->radialDiv(ii); + + double phiDelta = m_mgr->phiDelta(ii); + + bool fullPhiSector = false; + if (phiDelta == 0 || phiDelta >= 359.9 * CLHEP::degree) { + phiDelta = 360 * CLHEP::degree; + fullPhiSector = true; } + //else { + //phiDelta -= 2*phiepsilon; + //phiStart += phiepsilon; + // } + + if (shapeType == "UNKNOWN") { + if (radialDiv > 0) { + shapeType = "RADIAL"; + } else if (param->rmin() == rmin2 && param->rmax() == rmax2) { + if (fullPhiSector) { + shapeType = "TUBE"; + } else { + shapeType = "TUBS"; + } + } else { + shapeType = "CONS"; + } + } + + + int repeat = m_mgr->repeat(ii); + if (repeat == 0) repeat = 1; + + double phiStart = m_mgr->phiStart(ii); + double phiWidth = phiDelta; + + if (shapeType == "CONS" || shapeType == "TUBS") { + const double phiepsilon = 0.001 * CLHEP::degree; + phiWidth -= 2 * phiepsilon; + phiStart += phiepsilon; + } + + // Can be in degree or CLHEP::mm. Usually it is CLHEP::deg expect for BOX, TRAP and ROD shape + // Geometry manager makes no assumptions about units. So we must interpret here. + if (shapeType == "BOX" || shapeType == "ROD" || shapeType == "ROD2" || shapeType == "TRAP") { + phiWidth = m_mgr->width(ii); // in mm + } + + if (shapeType == "PGON" || shapeType == "PGON2" || + shapeType == "CONE" || shapeType == "CONS" || + shapeType == "PGON3" || shapeType == "PGON4") { + if ((rmin2 != param->rmin()) || (rmax2 != param->rmax())) { + needsRotation = true; + } + } + + int sides = 0; + int nCopies = 1; + if (shapeType == "PGON" || shapeType == "PGON2" || + shapeType == "PGON3" || shapeType == "PGON4") { + sides = repeat; + } else { + nCopies = repeat; + } + + // Force nCopies to 1 for TUBE and CONE + if (shapeType.empty() || shapeType == "TUBE" || shapeType == "CONE") { + nCopies = 1; + } + + param->setShapeType(shapeType); + param->setRmin2(rmin2); + param->setRmax2(rmax2); + param->setPhiLoc(phiStart); + param->setPhiWidth(phiWidth); + param->setSides(sides); + param->setNCopies(nCopies); + //param->setRadialDiv(radialDiv); + //param->setPhiStep(phiStep); } - - int sides = 0; - int nCopies = 1; - if (shapeType == "PGON" || shapeType == "PGON2" || - shapeType == "PGON3" || shapeType == "PGON4") { - sides = repeat; - } else { - nCopies = repeat; - } - - // Force nCopies to 1 for TUBE and CONE - if (shapeType.empty() || shapeType == "TUBE" || shapeType == "CONE") { - nCopies = 1; + + param->setNeedsRotation(needsRotation); + + + // + // If zin is 0... (within 10^-5) this is a volume symmetric around + // the origin + // + if (std::abs(param->zmin()) < 0.000001) { + param->setZmin(-param->zmax()); + param->setZsymm(false); } - - param->setShapeType(shapeType); - param->setRmin2(rmin2); - param->setRmax2(rmax2); - param->setPhiLoc(phiStart); - param->setPhiWidth(phiWidth); - param->setSides(sides); - param->setNCopies(nCopies); - //param->setRadialDiv(radialDiv); - //param->setPhiStep(phiStep); - } - - param->setNeedsRotation(needsRotation); - - - // - // If zin is 0... (within 10^-5) this is a volume symmetric around - // the origin - // - if(std::abs(param->zmin()) < 0.000001) { - param->setZmin(-param->zmax()); - param->setZsymm(false); - } - - param->setLabel(m_label, volId); - - return param; -} + param->setLabel(m_label, volId); + + return param; + } } // end namespace diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/TopLevelPlacements.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/TopLevelPlacements.cxx index 19645528405ab34be332d9156687a8f5752c7214..bd6b6f12ce86117c50ba3d599715e36c6b9da583 100755 --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/TopLevelPlacements.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/TopLevelPlacements.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/TopLevelPlacements.h" @@ -17,88 +17,80 @@ HepGeom::Transform3D TopLevelPlacements::s_identityTransform = HepGeom::Transfor TopLevelPlacements::TopLevelPlacements(IRDBRecordset_ptr topLevelTable) - : m_noTopLevelTable(true) -{ + : m_noTopLevelTable(true) { fillPlacements(topLevelTable); } -TopLevelPlacements::~TopLevelPlacements() -{ - std::map<std::string, Part *>::const_iterator iter; - for (iter=m_parts.begin(); iter != m_parts.end(); ++iter) delete iter->second; +TopLevelPlacements::~TopLevelPlacements() { + std::map<std::string, Part*>::const_iterator iter; + for (iter = m_parts.begin(); iter != m_parts.end(); ++iter) delete iter->second; } -const HepGeom::Transform3D & -TopLevelPlacements::transform(const std::string & partName) const -{ - Part * part = getPart(partName); +const HepGeom::Transform3D& +TopLevelPlacements::transform(const std::string& partName) const { + Part* part = getPart(partName); + if (part) return part->transform; + return s_identityTransform; } - -bool -TopLevelPlacements::present(const std::string & partName) const -{ +bool +TopLevelPlacements::present(const std::string& partName) const { // If no table present assume everything is present. if (m_noTopLevelTable) return true; - return (getPart(partName) != 0); + + return(getPart(partName) != 0); } void -TopLevelPlacements::fillPlacements(IRDBRecordset_ptr topLevelTable) -{ +TopLevelPlacements::fillPlacements(IRDBRecordset_ptr topLevelTable) { if (topLevelTable.get() == 0) { - m_noTopLevelTable = true; - return; + m_noTopLevelTable = true; + return; } m_noTopLevelTable = false; int numParts = topLevelTable->size(); for (int i = 0; i < numParts; i++) { - const IRDBRecord * record = (*topLevelTable)[i]; - std::string label = record->getString("LABEL"); + const IRDBRecord* record = (*topLevelTable)[i]; + std::string label = record->getString("LABEL"); - Part * part = new Part; + Part* part = new Part; part->label = label; part->transform = partTransform(record); - + m_parts[label] = part; - } - } - -HepGeom::Transform3D -TopLevelPlacements::partTransform(const IRDBRecord* record) const -{ - - double posX = record->getDouble("POSX") * CLHEP::mm; - double posY = record->getDouble("POSY") * CLHEP::mm; - double posZ = record->getDouble("POSZ") * CLHEP::mm; - double rotX = record->getDouble("ROTX") * CLHEP::degree; - double rotY = record->getDouble("ROTY") * CLHEP::degree; - double rotZ = record->getDouble("ROTZ") * CLHEP::degree; +HepGeom::Transform3D +TopLevelPlacements::partTransform(const IRDBRecord* record) const { + double posX = record->getDouble("POSX") * CLHEP::mm; + double posY = record->getDouble("POSY") * CLHEP::mm; + double posZ = record->getDouble("POSZ") * CLHEP::mm; + double rotX = record->getDouble("ROTX") * CLHEP::degree; + double rotY = record->getDouble("ROTY") * CLHEP::degree; + double rotZ = record->getDouble("ROTZ") * CLHEP::degree; int rotOrder = record->getInt("ROTORDER"); // Translation part - HepGeom::Transform3D transform = HepGeom::Translate3D(posX,posY,posZ); + HepGeom::Transform3D transform = HepGeom::Translate3D(posX, posY, posZ); // If rotation is zero return translation if (rotX == 0 && rotY == 0 && rotZ == 0) { return transform; } - // For rotation have to look at order. + // For rotation have to look at order. // 123 means rotate around X, then Y , then Z. // 312 means rotate around Z, then X , then Y. // etc - int ixyz1 = rotOrder/100 - 1; - int ixyz2 = (rotOrder%100) / 10 - 1; - int ixyz3 = (rotOrder%10) - 1 ; + int ixyz1 = rotOrder / 100 - 1; + int ixyz2 = (rotOrder % 100) / 10 - 1; + int ixyz3 = (rotOrder % 10) - 1; - if (ixyz1 < 0 || ixyz1 > 2 || + if (ixyz1 < 0 || ixyz1 > 2 || ixyz2 < 0 || ixyz2 > 2 || ixyz3 < 0 || ixyz3 > 2) { std::cout << "ERROR: Invalid rotation order:" << rotOrder << std::endl; @@ -106,9 +98,11 @@ TopLevelPlacements::partTransform(const IRDBRecord* record) const ixyz2 = 1; ixyz3 = 2; } - + // List of the three transforms - HepGeom::Transform3D * xformList[3] = {0, 0, 0}; + HepGeom::Transform3D* xformList[3] = { + 0, 0, 0 + }; if (rotX != 0) xformList[0] = new HepGeom::RotateX3D(rotX); if (rotY != 0) xformList[1] = new HepGeom::RotateY3D(rotY); if (rotZ != 0) xformList[2] = new HepGeom::RotateZ3D(rotZ); @@ -123,17 +117,13 @@ TopLevelPlacements::partTransform(const IRDBRecord* record) const delete xformList[2]; return transform * rotation; - } - -TopLevelPlacements::Part * -TopLevelPlacements::getPart(const std::string & partName) const -{ - std::map<std::string, Part *>::const_iterator iter; +TopLevelPlacements::Part* +TopLevelPlacements::getPart(const std::string& partName) const { + std::map<std::string, Part*>::const_iterator iter; iter = m_parts.find(partName); if (iter == m_parts.end()) return 0; + return iter->second; } - - diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/TubeVolData.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/TubeVolData.cxx index b5467756dc0bfb5961f9f683958892260a41741b..d12523d06b71de476a15857d63f657627bdc0341 100755 --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/TubeVolData.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/TubeVolData.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/TubeVolData.h" #include "RDBAccessSvc/IRDBRecord.h" @@ -11,25 +11,18 @@ #include <iostream> namespace InDetDD { - - - -std::string -TubeVolData::material() const -{ - return m_record->getString("MATERIAL"); -} - - -double -TubeVolData::maxRadius() const -{ - return std::max(m_rmax1,m_rmax2); -} - - -TubeVolData::TubeVolData(const IRDBRecord * record) - : m_record(record), + std::string + TubeVolData::material() const { + return m_record->getString("MATERIAL"); + } + + double + TubeVolData::maxRadius() const { + return std::max(m_rmax1, m_rmax2); + } + + TubeVolData::TubeVolData(const IRDBRecord* record) + : m_record(record), m_bothZ(false), m_nRepeat(0), m_radialDiv(0), @@ -41,75 +34,68 @@ TubeVolData::TubeVolData(const IRDBRecord * record) m_rmax1(0.), m_rmax2(0.), m_length(0.), - m_zMid(0.) -{ - // add an 2*epsilon gap between phi sectors. - const double phiepsilon = 0.001*CLHEP::degree; - - bool fullPhiSector = false; - - - // Get the parameters which we need to do some preprocessing with. - // The rest are obtained directly from RDB. + m_zMid(0.) { + // add an 2*epsilon gap between phi sectors. + const double phiepsilon = 0.001 * CLHEP::degree; + + bool fullPhiSector = false; + + + // Get the parameters which we need to do some preprocessing with. + // The rest are obtained directly from RDB. + + if (m_record) { + m_phiStart = m_record->getDouble("PHISTART") * CLHEP::degree; + m_phiDelta = m_record->getDouble("PHIDELTA") * CLHEP::degree; + m_phiStep = m_record->getDouble("PHISTEP") * CLHEP::degree; + m_nRepeat = m_record->getInt("NREPEAT"); + m_rmin1 = m_record->getDouble("RMIN") * CLHEP::mm; + m_rmax1 = m_record->getDouble("RMAX") * CLHEP::mm; + m_rmin2 = m_record->getDouble("RMIN2") * CLHEP::mm; + m_rmax2 = m_record->getDouble("RMAX2") * CLHEP::mm; + m_radialDiv = 0; + if (!m_record->isFieldNull("RADIAL")) { + m_radialDiv = m_record->getInt("RADIAL"); + } + m_bothZ = 0; + if (!m_record->isFieldNull("ZSYMM")) { + m_bothZ = m_record->getInt("ZSYMM"); + } - if(m_record) { - m_phiStart = m_record->getDouble("PHISTART") * CLHEP::degree; - m_phiDelta = m_record->getDouble("PHIDELTA") * CLHEP::degree; - m_phiStep = m_record->getDouble("PHISTEP") * CLHEP::degree; - m_nRepeat = m_record->getInt("NREPEAT"); - m_rmin1 = m_record->getDouble("RMIN") * CLHEP::mm; - m_rmax1 = m_record->getDouble("RMAX") * CLHEP::mm; - m_rmin2 = m_record->getDouble("RMIN2") * CLHEP::mm; - m_rmax2 = m_record->getDouble("RMAX2") * CLHEP::mm; - m_radialDiv = 0; - if (!m_record->isFieldNull("RADIAL")) { - m_radialDiv = m_record->getInt("RADIAL"); - } - m_bothZ = 0; - if (!m_record->isFieldNull("ZSYMM")) { - m_bothZ = m_record->getInt("ZSYMM"); - } + double zmin = m_record->getDouble("ZMIN") * CLHEP::mm; + double zmax = m_record->getDouble("ZMAX") * CLHEP::mm; + m_length = std::abs(zmax - zmin); + m_zMid = 0.5 * (zmin + zmax); - double zmin = m_record->getDouble("ZMIN") * CLHEP::mm; - double zmax = m_record->getDouble("ZMAX") * CLHEP::mm; - m_length = std::abs(zmax - zmin); - m_zMid = 0.5*(zmin+zmax); - - if (m_phiDelta == 0 || m_phiDelta >=359.9*CLHEP::degree) { - m_phiDelta = 360*CLHEP::degree; - fullPhiSector = true; - } else { - m_phiDelta -= 2*phiepsilon; - m_phiStart += phiepsilon; - } + if (m_phiDelta == 0 || m_phiDelta >= 359.9 * CLHEP::degree) { + m_phiDelta = 360 * CLHEP::degree; + fullPhiSector = true; + } else { + m_phiDelta -= 2 * phiepsilon; + m_phiStart += phiepsilon; + } - // Force nRepeat to be >= 1; - if (m_nRepeat <= 0) m_nRepeat = 1; - // if PHISTEP==0 then set it to be equi-distant steps filling up phi. - if (m_phiStep == 0) { - m_phiStep = 360*CLHEP::degree/m_nRepeat; - } + // Force nRepeat to be >= 1; + if (m_nRepeat <= 0) m_nRepeat = 1; + // if PHISTEP==0 then set it to be equi-distant steps filling up phi. + if (m_phiStep == 0) { + m_phiStep = 360 * CLHEP::degree / m_nRepeat; + } - if (m_rmin2 <= 0) m_rmin2 = m_rmin1; - if (m_rmax2 <= 0) m_rmax2 = m_rmax1; - - if (m_radialDiv > 0) { - m_shape = TubeVolData::RADIAL; - } else if (m_rmin1 == m_rmin2 && m_rmax1 == m_rmax2 ) { - if (fullPhiSector) { - m_shape = TubeVolData::TUBE; + if (m_rmin2 <= 0) m_rmin2 = m_rmin1; + if (m_rmax2 <= 0) m_rmax2 = m_rmax1; + + if (m_radialDiv > 0) { + m_shape = TubeVolData::RADIAL; + } else if (m_rmin1 == m_rmin2 && m_rmax1 == m_rmax2) { + if (fullPhiSector) { + m_shape = TubeVolData::TUBE; + } else { + m_shape = TubeVolData::TUBS; + } } else { - m_shape = TubeVolData::TUBS; + m_shape = TubeVolData::CONS; } - } else { - m_shape = TubeVolData::CONS; - } - } else std::cout << "Unexpected ERROR in ExtraMaterial!" << std::endl; - -} - + } else std::cout << "Unexpected ERROR in ExtraMaterial!" << std::endl; + } } // end namespace - - - - diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeBuilder.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeBuilder.cxx index bdfb6629ac3cca951a467d2385b75718c437ace6..8397564f503a6cf085464abb49799a48b025634a 100755 --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeBuilder.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeBuilder.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/VolumeBuilder.h" #include "InDetGeoModelUtils/ServiceVolume.h" @@ -10,445 +10,382 @@ #include "GeoModelKernel/GeoFullPhysVol.h" #include "GeoModelKernel/GeoMaterial.h" #include "GeoModelKernel/GeoTransform.h" - #include "CLHEP/Geometry/Transform3D.h" -//#include <sstream> -//#include <iomanip> -namespace InDetDD { -VolumeBuilder::VolumeBuilder(const Zone & zone, const std::vector<const ServiceVolume * > & services) - : m_msg("InDetDDVolumeBuilder"), +namespace InDetDD { + VolumeBuilder::VolumeBuilder(const Zone& zone, const std::vector<const ServiceVolume* >& services) + : m_msg("InDetDDVolumeBuilder"), m_region("None"), // Empty refers to a valid region. Set some default so we can check it is actually set. m_zcenter(0), m_services(0), m_servEnvelope(0), m_servChild(0), - m_matManager(0) -{ - m_splitter.splitAll(zone, services); -} + m_matManager(0) { + m_splitter.splitAll(zone, services); + } -VolumeBuilder::VolumeBuilder(const std::vector<const ServiceVolume * > & services) - : m_msg("InDetDDVolumeBuilder"), + VolumeBuilder::VolumeBuilder(const std::vector<const ServiceVolume* >& services) + : m_msg("InDetDDVolumeBuilder"), m_region("None"), // Empty refers to a valid region. Set some default so we can check it is actually set. m_zcenter(0), m_services(&services), m_servEnvelope(0), m_servChild(0), - - m_matManager(0) -{} - -VolumeBuilder::VolumeBuilder(const Zone & zone, const std::vector<const ServiceVolume * > & services, - const std::vector<const ServiceVolume * > & servEnvelope, - const std::vector<const ServiceVolume * > & servChild ) - : m_msg("InDetDDVolumeBuilder"), + {} + + VolumeBuilder::VolumeBuilder(const Zone& zone, const std::vector<const ServiceVolume* >& services, + const std::vector<const ServiceVolume* >& servEnvelope, + const std::vector<const ServiceVolume* >& servChild) + : m_msg("InDetDDVolumeBuilder"), m_region("None"), // Empty refers to a valid region. Set some default so we can check it is actually set. m_zcenter(0), m_services(0), m_servEnvelope(&servEnvelope), m_servChild(&servChild), - m_matManager(0) -{ - m_splitter.splitAll(zone, services); -} - -void VolumeBuilder::addServices(const Zone & zone, const std::vector<const ServiceVolume * > & services) -{ - m_splitter.splitAll(zone, services); -} - -void VolumeBuilder::setRegion(const std::string & region, double zcenter) -{ - m_region = region; - m_zcenter = zcenter; -} + m_matManager(0) { + m_splitter.splitAll(zone, services); + } -const std::vector<const ServiceVolume * > & VolumeBuilder::services() -{ - // Return volumes defined in VolumeSplitter - if (m_services) return *m_services; - return m_splitter.getVolumes(); -} + void + VolumeBuilder::addServices(const Zone& zone, const std::vector<const ServiceVolume* >& services) { + m_splitter.splitAll(zone, services); + } + + void + VolumeBuilder::setRegion(const std::string& region, double zcenter) { + m_region = region; + m_zcenter = zcenter; + } + + const std::vector<const ServiceVolume* >& + VolumeBuilder::services() { + // Return volumes defined in VolumeSplitter + if (m_services) return *m_services; + return m_splitter.getVolumes(); + } -const std::vector<const ServiceVolume * > & VolumeBuilder::servicesEnv() -{ - return *m_servEnvelope; -} + const std::vector<const ServiceVolume* >& + VolumeBuilder::servicesEnv() { + return *m_servEnvelope; + } -const std::vector<const ServiceVolume * > & VolumeBuilder::servicesChild() -{ - return *m_servChild; -} + const std::vector<const ServiceVolume* >& + VolumeBuilder::servicesChild() { + return *m_servChild; + } -void -VolumeBuilder::buildAndPlace(const std::string & region, GeoPhysVol * parent, double zcenter) { - // Get volumes defined by Volume splitter and add them on top GeoPhysVol - setRegion(region, zcenter); - for (unsigned int iElement = 0; iElement < services().size(); ++iElement) - if(!isEnvelopeOrChild(iElement)){ - GeoVPhysVol* physVol = build(iElement); - - if (physVol) { - physVol->ref(); - for (int iCopy = 0; iCopy < numCopies(iElement); ++iCopy) { - parent->add(getPlacement(iElement,iCopy)); - parent->add(physVol); + void + VolumeBuilder::buildAndPlace(const std::string& region, GeoPhysVol* parent, double zcenter) { + // Get volumes defined by Volume splitter and add them on top GeoPhysVol + setRegion(region, zcenter); + for (unsigned int iElement = 0; iElement < services().size(); ++iElement) + if (!isEnvelopeOrChild(iElement)) { + GeoVPhysVol* physVol = build(iElement); + if (physVol) { + physVol->ref(); + for (int iCopy = 0; iCopy < numCopies(iElement); ++iCopy) { + parent->add(getPlacement(iElement, iCopy)); + parent->add(physVol); + } + physVol->unref();//should delete if never added } - physVol->unref();//should delete if never added } - } -} + } - -void -VolumeBuilder::buildAndPlace(const std::string& region, GeoFullPhysVol* parent, double zcenter) { - // Get volumes defined by Volume splitter and add them on top GeoPhysVol - setRegion(region, zcenter); - for (unsigned int iElement = 0; iElement < services().size(); ++iElement) { - if(!isEnvelopeOrChild(iElement)){ - // GeoVPhysVol* physVol = build(iElement); - GeoVPhysVol* physVol = build(iElement); - if (physVol) { - physVol->ref(); - for (int iCopy = 0; iCopy < numCopies(iElement); ++iCopy) { - parent->add(getPlacement(iElement,iCopy)); - parent->add(physVol); + void + VolumeBuilder::buildAndPlace(const std::string& region, GeoFullPhysVol* parent, double zcenter) { + // Get volumes defined by Volume splitter and add them on top GeoPhysVol + setRegion(region, zcenter); + for (unsigned int iElement = 0; iElement < services().size(); ++iElement) { + if (!isEnvelopeOrChild(iElement)) { + // GeoVPhysVol* physVol = build(iElement); + GeoVPhysVol* physVol = build(iElement); + if (physVol) { + physVol->ref(); + for (int iCopy = 0; iCopy < numCopies(iElement); ++iCopy) { + parent->add(getPlacement(iElement, iCopy)); + parent->add(physVol); + } + physVol->unref();//should delete if never used } - physVol->unref();//should delete if never used } } - } - // if region is not Pixel -> stop here - if(region.compare("Pixel")!=0) return; - for (unsigned int iElement = 0; iElement < services().size(); ++iElement){ - if(getEnvelopeNum(iElement)>0&&services()[iElement]->envelopeParent()==0){ - buildAndPlaceEnvelope(region, parent, -1, iElement, zcenter); + // if region is not Pixel -> stop here + if (region.compare("Pixel") != 0) return; + + for (unsigned int iElement = 0; iElement < services().size(); ++iElement) { + if (getEnvelopeNum(iElement) > 0 && services()[iElement]->envelopeParent() == 0) { + buildAndPlaceEnvelope(region, parent, -1, iElement, zcenter); + } } } -} + void + VolumeBuilder::buildAndPlaceEnvelope(const std::string& region, GeoFullPhysVol* parent, int iParent, int iElement, + double zcenter) { + GeoPhysVol* physVol = dynamic_cast<GeoPhysVol*>(build(iElement)); -void -VolumeBuilder::buildAndPlaceEnvelope(const std::string& region, GeoFullPhysVol* parent, int iParent, int iElement, double zcenter){ - GeoPhysVol* physVol = dynamic_cast<GeoPhysVol*>(build(iElement)); - if (physVol) { - for (unsigned int iChild = 0; iChild < services().size(); ++iChild) { - if(isChildService(iElement,iChild)&&services()[iChild]->envelopeNum()>0){ - // if volume is a child volume : build and place it - buildAndPlaceEnvelope(region, physVol, iElement, iChild, zcenter); + if (physVol) { + for (unsigned int iChild = 0; iChild < services().size(); ++iChild) { + if (isChildService(iElement, iChild) && services()[iChild]->envelopeNum() > 0) { + // if volume is a child volume : build and place it + buildAndPlaceEnvelope(region, physVol, iElement, iChild, zcenter); + } } - } - for (unsigned int iChild = 0; iChild < services().size(); ++iChild) { - if(isChildService(iElement,iChild)&&services()[iChild]->envelopeNum()==0){ - // if volume is not a child volume - GeoVPhysVol* physVol_child = build(iChild); - if (physVol_child) { - physVol_child->ref(); - for (int iCopy2 = 0; iCopy2 < numCopies(iChild); ++iCopy2) { - physVol->add(getPlacementEnvelope(iChild,iCopy2,iElement)); - physVol->add(physVol_child); + for (unsigned int iChild = 0; iChild < services().size(); ++iChild) { + if (isChildService(iElement, iChild) && services()[iChild]->envelopeNum() == 0) { + // if volume is not a child volume + GeoVPhysVol* physVol_child = build(iChild); + if (physVol_child) { + physVol_child->ref(); + for (int iCopy2 = 0; iCopy2 < numCopies(iChild); ++iCopy2) { + physVol->add(getPlacementEnvelope(iChild, iCopy2, iElement)); + physVol->add(physVol_child); + } + physVol_child->unref();//should delete if was never added } - physVol_child->unref();//should delete if was never added } } - } - - for (int iCopy = 0; iCopy < numCopies(iElement); ++iCopy) { - // add all the copies - if(iParent<0)parent->add(getPlacement(iElement,iCopy)); - else parent->add(getPlacementEnvelope(iElement,iCopy,iParent)); - parent->add(physVol); + for (int iCopy = 0; iCopy < numCopies(iElement); ++iCopy) { + // add all the copies + if (iParent < 0) parent->add(getPlacement(iElement, iCopy)); + else parent->add(getPlacementEnvelope(iElement, iCopy, iParent)); + parent->add(physVol); + } } } -} -void -VolumeBuilder::buildAndPlaceEnvelope(const std::string& region, GeoPhysVol* parent, int iParent, int iElement, double zcenter){ - GeoPhysVol* physVol = dynamic_cast<GeoPhysVol*>(build(iElement)); - if (physVol){ - physVol->ref(); - for (unsigned int iChild = 0; iChild < services().size(); ++iChild) { - if(isChildService(iElement,iChild)&&services()[iChild]->envelopeNum()>0){ - // if volume is a child volume : build and place it - buildAndPlaceEnvelope(region, physVol, iElement, iChild, zcenter); + void + VolumeBuilder::buildAndPlaceEnvelope(const std::string& region, GeoPhysVol* parent, int iParent, int iElement, + double zcenter) { + GeoPhysVol* physVol = dynamic_cast<GeoPhysVol*>(build(iElement)); + if (physVol) { + physVol->ref(); + for (unsigned int iChild = 0; iChild < services().size(); ++iChild) { + if (isChildService(iElement, iChild) && services()[iChild]->envelopeNum() > 0) { + // if volume is a child volume : build and place it + buildAndPlaceEnvelope(region, physVol, iElement, iChild, zcenter); + } } - } - for (unsigned int iChild = 0; iChild < services().size(); ++iChild) { - if(isChildService(iElement,iChild)&&services()[iChild]->envelopeNum()==0){ - // if volume is not a child volume - GeoVPhysVol* physVol_child = build(iChild); - if (physVol_child) { - physVol_child->ref(); - for (int iCopy2 = 0; iCopy2 < numCopies(iChild); ++iCopy2) { - physVol->add(getPlacementEnvelope(iChild,iCopy2,iElement)); - physVol->add(physVol_child); + for (unsigned int iChild = 0; iChild < services().size(); ++iChild) { + if (isChildService(iElement, iChild) && services()[iChild]->envelopeNum() == 0) { + // if volume is not a child volume + GeoVPhysVol* physVol_child = build(iChild); + if (physVol_child) { + physVol_child->ref(); + for (int iCopy2 = 0; iCopy2 < numCopies(iChild); ++iCopy2) { + physVol->add(getPlacementEnvelope(iChild, iCopy2, iElement)); + physVol->add(physVol_child); + } + physVol_child->unref(); //will delete physVol_child if never used } - physVol_child->unref(); //will delete physVol_child if never used } } + for (int iCopy = 0; iCopy < numCopies(iElement); ++iCopy) { + // add all the copies + if (iParent < 0) parent->add(getPlacement(iElement, iCopy)); + else parent->add(getPlacementEnvelope(iElement, iCopy, iParent)); + parent->add(physVol); + } + physVol->unref();//will delete physvol if it was never used } - for (int iCopy = 0; iCopy < numCopies(iElement); ++iCopy) { - // add all the copies - if(iParent<0)parent->add(getPlacement(iElement,iCopy)); - else parent->add(getPlacementEnvelope(iElement,iCopy,iParent)); - parent->add(physVol); - } - physVol->unref();//will delete physvol if it was never used } -} -GeoVPhysVol* -VolumeBuilder::build(int iElement) { - if (m_region == "None") { - msg(MSG::ERROR) << "No region set. Cannot build services"<< endmsg; - return nullptr; - } - const ServiceVolume & param = *(services()[iElement]); - // If the subelement does not belong to the current region return 0. - if (param.region() != m_region) return 0; - const GeoShape * serviceShape = param.getShape(); - double volume = param.origVolume(); - std::string logName = param.fullLabel(); - const GeoMaterial* serviceMat = param.material(); - std::string materialName; - if (!serviceMat) { - materialName = param.materialName(); - if (m_matManager){ - //serviceMat = m_matManager->getMaterialForVolume(materialName,volume/param.fractionInRegion()); - // FIXME - serviceMat = m_matManager->getMaterialForVolume(materialName,volume); + GeoVPhysVol* + VolumeBuilder::build(int iElement) { + if (m_region == "None") { + msg(MSG::ERROR) << "No region set. Cannot build services" << endmsg; + return nullptr; + } + const ServiceVolume& param = *(services()[iElement]); + // If the subelement does not belong to the current region return 0. + if (param.region() != m_region) return 0; + + const GeoShape* serviceShape = param.getShape(); + double volume = param.origVolume(); + std::string logName = param.fullLabel(); + const GeoMaterial* serviceMat = param.material(); + std::string materialName; + if (!serviceMat) { + materialName = param.materialName(); + if (m_matManager) { + //serviceMat = m_matManager->getMaterialForVolume(materialName,volume/param.fractionInRegion()); + // FIXME + serviceMat = m_matManager->getMaterialForVolume(materialName, volume); + } else { + msg(MSG::ERROR) << "Material manager not available. Cannot build material." << endmsg; + return nullptr; + } } else { - msg(MSG::ERROR) << "Material manager not available. Cannot build material." << endmsg; - return nullptr ; - } - } else { - materialName = serviceMat->getName(); - } - if (msgLvl(MSG::DEBUG)) { - msg(MSG::DEBUG) << "Volume/material: " << logName << "/" << materialName << endmsg; - if (!param.shapeType().empty()) msg(MSG::DEBUG) << " shape: " << param.shapeType() << endmsg; - msg(MSG::DEBUG) << " volume (CLHEP::cm3): " << volume/CLHEP::cm3 << endmsg; - msg(MSG::DEBUG) << " rmin,rmax,zmin,zmax: " - << param.rmin() << ", " - << param.rmax() << ", " - << param.zmin() << ", " - << param.zmax() << endmsg; + materialName = serviceMat->getName(); + } + if (msgLvl(MSG::DEBUG)) { + msg(MSG::DEBUG) << "Volume/material: " << logName << "/" << materialName << endmsg; + if (!param.shapeType().empty()) msg(MSG::DEBUG) << " shape: " << param.shapeType() << endmsg; + msg(MSG::DEBUG) << " volume (CLHEP::cm3): " << volume / CLHEP::cm3 << endmsg; + msg(MSG::DEBUG) << " rmin,rmax,zmin,zmax: " + << param.rmin() << ", " + << param.rmax() << ", " + << param.zmin() << ", " + << param.zmax() << endmsg; + } + // Or use volume of original volume in param. + //const GeoMaterial* serviceMat = mat_mgr->getMaterialForVolume(param.material(),param.origVolume()); + GeoLogVol* serviceLog = new GeoLogVol(logName, serviceShape, serviceMat); + GeoPhysVol* servicePhys = new GeoPhysVol(serviceLog); + return servicePhys; } - // Or use volume of original volume in param. - //const GeoMaterial* serviceMat = mat_mgr->getMaterialForVolume(param.material(),param.origVolume()); - GeoLogVol* serviceLog = new GeoLogVol(logName,serviceShape,serviceMat); - GeoPhysVol* servicePhys = new GeoPhysVol(serviceLog); - return servicePhys; -} - - -bool VolumeBuilder::isEnvelopeOrChild(int iElement){ - - const ServiceVolume& param = *(services()[iElement]); - if(param.envelopeNum()==0&¶m.envelopeParent()==0)return false; - return true; -} -int VolumeBuilder::getEnvelopeNum(int iElement){ - const ServiceVolume& param = *(services()[iElement]); - return param.envelopeNum(); -} + bool + VolumeBuilder::isEnvelopeOrChild(int iElement) { + const ServiceVolume& param = *(services()[iElement]); + if (param.envelopeNum() == 0 && param.envelopeParent() == 0) return false; + return true; + } -int VolumeBuilder::getParentNum(int iElement){ - const ServiceVolume& param = *(services()[iElement]); - return param.envelopeParent(); -} + int + VolumeBuilder::getEnvelopeNum(int iElement) { + const ServiceVolume& param = *(services()[iElement]); + return param.envelopeNum(); + } -double VolumeBuilder::getZcenter(int iElement){ - const ServiceVolume& param = *(services()[iElement]); - return (param.zmin()+param.zmax())*0.5; -} + int + VolumeBuilder::getParentNum(int iElement) { + const ServiceVolume& param = *(services()[iElement]); + return param.envelopeParent(); + } -bool VolumeBuilder::isChildService(int iElt,int iChld) -{ - const ServiceVolume& param1 = *(services()[iElt]); - const ServiceVolume& param2 = *(services()[iChld]); - - if(iElt==iChld||param1.envelopeNum()!=param2.envelopeParent())return false; + double + VolumeBuilder::getZcenter(int iElement) { + const ServiceVolume& param = *(services()[iElement]); + return (param.zmin() + param.zmax()) * 0.5; + } - if(param1.zsymm()==1) - { - double zmin=(param1.zmin()*param2.zmin()); - double zmax=(param1.zmax()*param2.zmax()); - - if(zmin>0&&zmax>0) return true; + bool + VolumeBuilder::isChildService(int iElt, int iChld) { + const ServiceVolume& param1 = *(services()[iElt]); + const ServiceVolume& param2 = *(services()[iChld]); + if (iElt == iChld || param1.envelopeNum() != param2.envelopeParent()) return false; + if (param1.zsymm() == 1) { + double zmin = (param1.zmin() * param2.zmin()); + double zmax = (param1.zmax() * param2.zmax()); + if (zmin > 0 && zmax > 0) return true; return false; } + return true; + } -// double zmin=(param1.zmin()*param2.zmin()); -// double zmax=(param1.zmax()*param2.zmax()); - -// if(zmin>0&&zmax>0) return true; - return true; - -} - -int VolumeBuilder::numCopies(int iElement) -{ - return services()[iElement]->nCopies(); -} - - -GeoTransform * VolumeBuilder::getPlacement(int iElement, int iCopy) -{ - const ServiceVolume & param = *(services()[iElement]); - - // NB. Corrected for placement in endcaps - double zpos = param.zposition() - m_zcenter; - - // Shift along Z axis ( IBL shift ) - double zshift = param.zShift(); - zpos += zshift; - - // Check if we need to rotate around Y axis. - bool rotateAroundY = false; - if (param.needsRotation()) { // zpos will always be negative in this case - zpos = -zpos; - rotateAroundY = true; + int + VolumeBuilder::numCopies(int iElement) { + return services()[iElement]->nCopies(); } - - // Most are just translated in z - HepGeom::Transform3D xform = HepGeom::TranslateZ3D(zpos); - - double phiStart = 0; - - // BOX, ROD and TRAP need special treatment. - const std::string & shapeType = param.shapeType(); - if (shapeType == "TRAP" || shapeType == "TRAP2") - { - // Need to rotate by -90 deg. - xform = HepGeom::RotateZ3D(-90.*CLHEP::deg) * xform; + + GeoTransform* + VolumeBuilder::getPlacement(int iElement, int iCopy) { + const ServiceVolume& param = *(services()[iElement]); + // NB. Corrected for placement in endcaps + double zpos = param.zposition() - m_zcenter; + // Shift along Z axis ( IBL shift ) + double zshift = param.zShift(); + zpos += zshift; + // Check if we need to rotate around Y axis. + bool rotateAroundY = false; + if (param.needsRotation()) { // zpos will always be negative in this case + zpos = -zpos; + rotateAroundY = true; } - if ( shapeType == "TRAP2") - { + // Most are just translated in z + HepGeom::Transform3D xform = HepGeom::TranslateZ3D(zpos); + double phiStart = 0; + // BOX, ROD and TRAP need special treatment. + const std::string& shapeType = param.shapeType(); + if (shapeType == "TRAP" || shapeType == "TRAP2") { // Need to rotate by -90 deg. - // xform = HepGeom::RotateZ3D(-90.*CLHEP::deg) * HepGeom::RotateY3D(-90.*CLHEP::deg) * xform; - - xform = HepGeom::RotateZ3D(-90.*CLHEP::deg) * xform; // * HepGeom::RotateX3D(-90.*CLHEP::deg); + xform = HepGeom::RotateZ3D(-90. * CLHEP::deg) * xform; } - if (shapeType == "BOX" || shapeType == "TRAP" || shapeType == "TRAP2") - { - double radius = 0.5*(param.rmin() + param.rmax()); + if (shapeType == "TRAP2") { + xform = HepGeom::RotateZ3D(-90. * CLHEP::deg) * xform; // * HepGeom::RotateX3D(-90.*CLHEP::deg); + } + if (shapeType == "BOX" || shapeType == "TRAP" || shapeType == "TRAP2") { + double radius = 0.5 * (param.rmin() + param.rmax()); xform = HepGeom::TranslateX3D(radius) * xform; - phiStart = param.phiLoc(); - } - else if (shapeType == "ROD" || shapeType == "ROD2") - { + phiStart = param.phiLoc(); + } else if (shapeType == "ROD" || shapeType == "ROD2") { double radius = param.rmin(); xform = HepGeom::TranslateX3D(radius) * xform; - phiStart = param.phiLoc(); + phiStart = param.phiLoc(); } - - // For volumes that are placed more than once. - double deltaPhi = 0; - if (iCopy > 0) { - deltaPhi = 2.*M_PI/param.nCopies(); - } - - double phi = phiStart + deltaPhi * iCopy; - if (phi) { - xform = HepGeom::RotateZ3D(phi) * xform; - } - - // For shapes that are not symmetric about a rotation around Y axis. We need to rotate. - if (rotateAroundY) { - xform = HepGeom::RotateY3D(180.*CLHEP::degree) * xform; - } - - return new GeoTransform(xform); -} - - -GeoTransform * VolumeBuilder::getPlacementEnvelope(int iElement, int iCopy, int iEnvElement) -{ - const ServiceVolume & param = *(services()[iElement]); - - const ServiceVolume& paramEnv = *(services()[iEnvElement]); - double zCenter=(paramEnv.zmin()+paramEnv.zmax())*0.5; - double rCenter=0.; - - bool bMoveToCenter=false; - if(paramEnv.shapeType()=="BOX")bMoveToCenter=true; - if(paramEnv.shapeType()=="TUBE"&¶mEnv.zsymm()==1&&fabs(paramEnv.zmin())>0.01)bMoveToCenter=true; - // if(paramEnv.shapeType()=="PGON31")bMoveToCenter=true; - if(bMoveToCenter)rCenter=(paramEnv.rmin()+paramEnv.rmax())*0.5; - - // std::cout<<"Placement env : "<<param.volName()<<" "<<paramEnv.volName()<<" "<<rCenter<<std::endl; - - // NB. Corrected for placement in endcaps - double zpos = param.zposition() - zCenter; - // double rpos = param.rposition() - rCenter; - - // Check if we need to rotate around Y axis. - bool rotateAroundY = false; - if (param.needsRotation()) { // zpos will always be negative in this case - zpos = -zpos; - rotateAroundY = true; + // For volumes that are placed more than once. + double deltaPhi = 0; + if (iCopy > 0) { + deltaPhi = 2. * M_PI / param.nCopies(); + } + double phi = phiStart + deltaPhi * iCopy; + if (phi) { + xform = HepGeom::RotateZ3D(phi) * xform; + } + // For shapes that are not symmetric about a rotation around Y axis. We need to rotate. + if (rotateAroundY) { + xform = HepGeom::RotateY3D(180. * CLHEP::degree) * xform; + } + return new GeoTransform(xform); } - // Most are just translated in z - HepGeom::Transform3D xform = HepGeom::TranslateZ3D(zpos); - - const std::string & shapeType = param.shapeType(); - // if (shapeType == "BOX") - // xform = HepGeom::TranslateY3D(rpos) *HepGeom::TranslateZ3D(zpos); - - double phiStart = 0; - - // BOX, ROD and TRAP need special treatment. - // const std::string & shapeType = param.shapeType(); - if (shapeType == "TRAP" ) - { + GeoTransform* + VolumeBuilder::getPlacementEnvelope(int iElement, int iCopy, int iEnvElement) { + const ServiceVolume& param = *(services()[iElement]); + const ServiceVolume& paramEnv = *(services()[iEnvElement]); + double zCenter = (paramEnv.zmin() + paramEnv.zmax()) * 0.5; + double rCenter = 0.; + bool bMoveToCenter = false; + if (paramEnv.shapeType() == "BOX") bMoveToCenter = true; + if (paramEnv.shapeType() == "TUBE" && paramEnv.zsymm() == 1 && fabs(paramEnv.zmin()) > 0.01) bMoveToCenter = true; + if (bMoveToCenter) rCenter = (paramEnv.rmin() + paramEnv.rmax()) * 0.5; + // NB. Corrected for placement in endcaps + double zpos = param.zposition() - zCenter; + // Check if we need to rotate around Y axis. + bool rotateAroundY = false; + if (param.needsRotation()) { // zpos will always be negative in this case + zpos = -zpos; + rotateAroundY = true; + } + // Most are just translated in z + HepGeom::Transform3D xform = HepGeom::TranslateZ3D(zpos); + const std::string& shapeType = param.shapeType(); + double phiStart = 0; + // BOX, ROD and TRAP need special treatment. + if (shapeType == "TRAP") { // Need to rotate by -90 deg. - xform = HepGeom::RotateZ3D(-90.*CLHEP::deg) * xform; + xform = HepGeom::RotateZ3D(-90. * CLHEP::deg) * xform; } - if ( shapeType == "TRAP2") - { + if (shapeType == "TRAP2") { // Need to rotate by -90 deg. - // xform = HepGeom::RotateZ3D(-90.*CLHEP::deg) * HepGeom::RotateY3D(-90.*CLHEP::deg) * xform; - - xform = HepGeom::RotateX3D(-90.*CLHEP::deg) * xform; + xform = HepGeom::RotateX3D(-90. * CLHEP::deg) * xform; } - if (shapeType == "BOX" || shapeType == "TRAP"|| shapeType == "TRAP2") - { - double radius = 0.5*(param.rmin() + param.rmax()) - rCenter; + if (shapeType == "BOX" || shapeType == "TRAP" || shapeType == "TRAP2") { + double radius = 0.5 * (param.rmin() + param.rmax()) - rCenter; xform = HepGeom::TranslateX3D(radius) * xform; - phiStart = param.phiLoc(); - } - else if (shapeType == "ROD" || shapeType == "ROD2") - { + phiStart = param.phiLoc(); + } else if (shapeType == "ROD" || shapeType == "ROD2") { double radius = param.rmin(); xform = HepGeom::TranslateX3D(radius) * xform; - phiStart = param.phiLoc(); + phiStart = param.phiLoc(); } - - // For volumes that are placed more than once. - double deltaPhi = 0; - if (iCopy > 0) { - deltaPhi = 2.*M_PI/param.nCopies(); - } - - double phi = phiStart + deltaPhi * iCopy; - if (phi) { - xform = HepGeom::RotateZ3D(phi) * xform; - } - - // For shapes that are not symmetric about a rotation around Y axis. We need to rotate. - if (rotateAroundY) { - xform = HepGeom::RotateY3D(180.*CLHEP::degree) * xform; + // For volumes that are placed more than once. + double deltaPhi = 0; + if (iCopy > 0) { + deltaPhi = 2. * M_PI / param.nCopies(); + } + double phi = phiStart + deltaPhi * iCopy; + if (phi) { + xform = HepGeom::RotateZ3D(phi) * xform; + } + // For shapes that are not symmetric about a rotation around Y axis. We need to rotate. + if (rotateAroundY) { + xform = HepGeom::RotateY3D(180. * CLHEP::degree) * xform; + } + return new GeoTransform(xform); } - - return new GeoTransform(xform); -} - } // end namespace - - diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeSplitter.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeSplitter.cxx old mode 100644 new mode 100755 index e0ef2f76faa460816b9c41e41fabe9fb5a99bc7f..35fb29ef43a971a7cf6a547a0c9f94ef59f56b1d --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeSplitter.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeSplitter.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/VolumeSplitter.h" #include "InDetGeoModelUtils/VolumeSplitterUtils.h" @@ -11,191 +11,181 @@ #include <algorithm> namespace InDetDD { - -VolumeSplitter::VolumeSplitter() - : m_ownVolumes(true), - m_epsilon(0.0001*CLHEP::mm) -{} - -VolumeSplitter::~VolumeSplitter() -{ - if (m_ownVolumes) { - for (std::vector<const ServiceVolume *>::iterator iter = m_volumes.begin(); iter != m_volumes.end(); ++iter) { - delete *iter; + VolumeSplitter::VolumeSplitter() + : m_ownVolumes(true), + m_epsilon(0.0001 * CLHEP::mm) + {} + + VolumeSplitter::~VolumeSplitter() { + if (m_ownVolumes) { + for (std::vector<const ServiceVolume*>::iterator iter = m_volumes.begin(); iter != m_volumes.end(); ++iter) { + delete *iter; + } } } -} -const std::vector<const ServiceVolume *> & -VolumeSplitter::splitAll(const Zone & zone, const std::vector<const ServiceVolume *> origVolumeList) -{ - for (unsigned int i = 0; i < origVolumeList.size(); ++i) { - split(zone, *(origVolumeList[i])); - } - return m_volumes; -} - -void -VolumeSplitter::split(const Zone & zone, const ServiceVolume & origVolume) -{ - makeVolumes(&zone, origVolume); -} - -Ray -VolumeSplitter::makeRay(const ServiceVolume & origVolume) -{ - double zmin = origVolume.zmin(); - double zmax = origVolume.zmax(); - double rmin = origVolume.rmin(); - double rmax = origVolume.rmax(); - - if (zmax < zmin) std::swap(zmin,zmax); - if (rmax < rmin) std::swap(rmin,rmax); - - bool horizontal = (std::abs(zmax - zmin) > std::abs(rmax - rmin)); - if (horizontal) { - double rmid = 0.5*(rmin+rmax); - return Ray(Point(zmin,rmid), Point(zmax,rmid)); - } else { - double zmid = 0.5*(zmin+zmax); - return Ray(Point(zmid,rmin), Point(zmid,rmax)); + const std::vector<const ServiceVolume*>& + VolumeSplitter::splitAll(const Zone& zone, const std::vector<const ServiceVolume*> origVolumeList) { + for (unsigned int i = 0; i < origVolumeList.size(); ++i) { + split(zone, *(origVolumeList[i])); + } + return m_volumes; } -} - - -void -VolumeSplitter::makeVolumes(const Zone * zone, const ServiceVolume & origVolume) -{ - // See if its a symmetric volume - if (origVolume.zsymm()) { - ServiceVolume param1(origVolume); - ServiceVolume param2(origVolume); - param1.addLabel("A"); - param1.setNeedsRotation(false); - param2.addLabel("C"); - param2.setZmin(-origVolume.zmax()); - param2.setZmax(-origVolume.zmin()); - param2.setNeedsRotation(origVolume.needsRotation()); - splitVolume(zone, param1); - splitVolume(zone, param2); - } else { - ServiceVolume param1(origVolume); - param1.setNeedsRotation(false); - splitVolume(zone, param1); + + void + VolumeSplitter::split(const Zone& zone, const ServiceVolume& origVolume) { + makeVolumes(&zone, origVolume); } -} - -void -VolumeSplitter::splitVolume(const Zone * zone, const ServiceVolume & param) -{ - Ray ray = makeRay(param); - SegmentSplitter splitter; - const SegmentList & segments = splitter.split(zone, ray); - - //std::cout << "Segments: " << std::endl; - //segments.print(); - - double volumeTotal = 0; - std::vector<ServiceVolume * > tmpServiceVec; - tmpServiceVec.reserve(segments.size()); - for (unsigned int i = 0; i < segments.size(); ++i) { - const Segment & seg = segments.getSegment(i); - ServiceVolume * paramNew = new ServiceVolume(param); - - if (segments.horizontal()) { - if (param.splittableInZ()) { - paramNew->setZmin(seg.zmin()); - paramNew->setZmax(seg.zmax()); - // If z changed adjust r in cone like volumes - adjustR(param, paramNew); - } else if (i == 1) { - std::cout << "Volume " << param.fullLabel() << " cannot be split in Z" << std::endl; - } + + Ray + VolumeSplitter::makeRay(const ServiceVolume& origVolume) { + double zmin = origVolume.zmin(); + double zmax = origVolume.zmax(); + double rmin = origVolume.rmin(); + double rmax = origVolume.rmax(); + + if (zmax < zmin) std::swap(zmin, zmax); + if (rmax < rmin) std::swap(rmin, rmax); + + bool horizontal = (std::abs(zmax - zmin) > std::abs(rmax - rmin)); + if (horizontal) { + double rmid = 0.5 * (rmin + rmax); + return Ray(Point(zmin, rmid), Point(zmax, rmid)); } else { - // vertical - if (param.splittableInR()) { - paramNew->setRmin(seg.rmin()); - paramNew->setRmax(seg.rmax()); - } else if (i == 1) { - std::cout << "Volume " << param.fullLabel() << " cannot be split in R" << std::endl; - } + double zmid = 0.5 * (zmin + zmax); + return Ray(Point(zmid, rmin), Point(zmid, rmax)); } + } - if (seg.rotated()) { - // Often the -ve endcap region is created as the +ve endcap and then rotated. Therefore we have to reflect the z coords. - // Also have to turn off the needsRotation flag in case it was set. - //std::cout << "Segment " << i << " rotated " << paramNew->zmin() << " " << paramNew->zmax() <<std::endl; - double zminTmp = paramNew->zmin(); - paramNew->setZmin(-paramNew->zmax()); - paramNew->setZmax(-zminTmp); - paramNew->setNeedsRotation(false); - //std::cout << "After adjusted " << paramNew->zmin() << " " << paramNew->zmax() <<std::endl; + void + VolumeSplitter::makeVolumes(const Zone* zone, const ServiceVolume& origVolume) { + // See if its a symmetric volume + if (origVolume.zsymm()) { + ServiceVolume param1(origVolume); + ServiceVolume param2(origVolume); + param1.addLabel("A"); + param1.setNeedsRotation(false); + param2.addLabel("C"); + param2.setZmin(-origVolume.zmax()); + param2.setZmax(-origVolume.zmin()); + param2.setNeedsRotation(origVolume.needsRotation()); + splitVolume(zone, param1); + splitVolume(zone, param2); + } else { + ServiceVolume param1(origVolume); + param1.setNeedsRotation(false); + splitVolume(zone, param1); } + } - if (param.splittableInZ()) { - paramNew->reduceSize(m_epsilon); // safety gap between volumes - } + void + VolumeSplitter::splitVolume(const Zone* zone, const ServiceVolume& param) { + Ray ray = makeRay(param); + SegmentSplitter splitter; + const SegmentList& segments = splitter.split(zone, ray); + + //std::cout << "Segments: " << std::endl; + //segments.print(); + + double volumeTotal = 0; + + std::vector<ServiceVolume* > tmpServiceVec; + tmpServiceVec.reserve(segments.size()); + for (unsigned int i = 0; i < segments.size(); ++i) { + const Segment& seg = segments.getSegment(i); + ServiceVolume* paramNew = new ServiceVolume(param); + + if (segments.horizontal()) { + if (param.splittableInZ()) { + paramNew->setZmin(seg.zmin()); + paramNew->setZmax(seg.zmax()); + // If z changed adjust r in cone like volumes + adjustR(param, paramNew); + } else if (i == 1) { + std::cout << "Volume " << param.fullLabel() << " cannot be split in Z" << std::endl; + } + } else { + // vertical + if (param.splittableInR()) { + paramNew->setRmin(seg.rmin()); + paramNew->setRmax(seg.rmax()); + } else if (i == 1) { + std::cout << "Volume " << param.fullLabel() << " cannot be split in R" << std::endl; + } + } - paramNew->setRegion(seg.label()); - tmpServiceVec.push_back(paramNew); - - paramNew->getShape(); - volumeTotal += paramNew->volume(); - //std::cout << i << ": volume: " << paramNew->volume() << std::endl; + if (seg.rotated()) { + // Often the -ve endcap region is created as the +ve endcap and then rotated. Therefore we have to reflect the z + // coords. + // Also have to turn off the needsRotation flag in case it was set. + //std::cout << "Segment " << i << " rotated " << paramNew->zmin() << " " << paramNew->zmax() <<std::endl; + double zminTmp = paramNew->zmin(); + paramNew->setZmin(-paramNew->zmax()); + paramNew->setZmax(-zminTmp); + paramNew->setNeedsRotation(false); + //std::cout << "After adjusted " << paramNew->zmin() << " " << paramNew->zmax() <<std::endl; + } - } - //std::cout << "Total volume: " << volumeTotal << std::endl; - for (unsigned int i = 0; i < segments.size(); ++i) { - ServiceVolume * paramNew = tmpServiceVec[i]; - if (segments.size() > 1) { - std::ostringstream ostr; - ostr << "_" << i+1; - paramNew->addLabel(ostr.str()); - } - paramNew->setOrigVolume(volumeTotal); - m_volumes.push_back(paramNew); - } + if (param.splittableInZ()) { + paramNew->reduceSize(m_epsilon); // safety gap between volumes + } -} + paramNew->setRegion(seg.label()); + tmpServiceVec.push_back(paramNew); -// This takes care of cone like volumes and adjust the CLHEP::radius according to the adjusted z. -void -VolumeSplitter::adjustR(const ServiceVolume & param, ServiceVolume * paramNew) -{ - double z1 = param.zmin(); - double z2 = param.zmax(); - double z1New = paramNew->zmin(); - double z2New = paramNew->zmax(); - double rmin1 = param.rmin(); - double rmin2 = param.rmin2(); - double rmax1 = param.rmax(); - double rmax2 = param.rmax2(); - - // Make sure z1 < z2 - // as r1, r2 order assumes this - if (z1 > z2) std::swap(z1,z2); - if (z1New > z2New) std::swap(z1New,z2New); - - // Avoid divide by zero. - if (z1 == z2) return; - - if (z1 != z1New || z2 != z2New) { - if (rmin2 > 0 && rmax2 > 0 && (rmin1 != rmin2 || rmax1 != rmax2)) { - double slopeMin = (rmin2-rmin1)/(z2-z1); - double slopeMax = (rmax2-rmax1)/(z2-z1); - - double rmin1New = (z1New - z1) * slopeMin + rmin1; - double rmin2New = (z2New - z1) * slopeMin + rmin1; - - double rmax1New = (z1New - z1) * slopeMax + rmax1; - double rmax2New = (z2New - z1) * slopeMax + rmax1; - - paramNew->setRmin(rmin1New); - paramNew->setRmax(rmax1New); - paramNew->setRmin2(rmin2New); - paramNew->setRmax2(rmax2New); + paramNew->getShape(); + volumeTotal += paramNew->volume(); + //std::cout << i << ": volume: " << paramNew->volume() << std::endl; + } + //std::cout << "Total volume: " << volumeTotal << std::endl; + for (unsigned int i = 0; i < segments.size(); ++i) { + ServiceVolume* paramNew = tmpServiceVec[i]; + if (segments.size() > 1) { + std::ostringstream ostr; + ostr << "_" << i + 1; + paramNew->addLabel(ostr.str()); + } + paramNew->setOrigVolume(volumeTotal); + m_volumes.push_back(paramNew); } } -} +// This takes care of cone like volumes and adjust the CLHEP::radius according to the adjusted z. + void + VolumeSplitter::adjustR(const ServiceVolume& param, ServiceVolume* paramNew) { + double z1 = param.zmin(); + double z2 = param.zmax(); + double z1New = paramNew->zmin(); + double z2New = paramNew->zmax(); + double rmin1 = param.rmin(); + double rmin2 = param.rmin2(); + double rmax1 = param.rmax(); + double rmax2 = param.rmax2(); + + // Make sure z1 < z2 + // as r1, r2 order assumes this + if (z1 > z2) std::swap(z1, z2); + if (z1New > z2New) std::swap(z1New, z2New); + + // Avoid divide by zero. + if (z1 == z2) return; + + if (z1 != z1New || z2 != z2New) { + if (rmin2 > 0 && rmax2 > 0 && (rmin1 != rmin2 || rmax1 != rmax2)) { + double slopeMin = (rmin2 - rmin1) / (z2 - z1); + double slopeMax = (rmax2 - rmax1) / (z2 - z1); + + double rmin1New = (z1New - z1) * slopeMin + rmin1; + double rmin2New = (z2New - z1) * slopeMin + rmin1; + + double rmax1New = (z1New - z1) * slopeMax + rmax1; + double rmax2New = (z2New - z1) * slopeMax + rmax1; + + paramNew->setRmin(rmin1New); + paramNew->setRmax(rmax1New); + paramNew->setRmin2(rmin2New); + paramNew->setRmax2(rmax2New); + } + } + } } // end namespace diff --git a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeSplitterUtils.cxx b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeSplitterUtils.cxx old mode 100644 new mode 100755 index 887932cc537d444aba8f94fe6a72416e7803065b..4e0c4173df819b0f1475c2a6f4aac64620d48a18 --- a/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeSplitterUtils.cxx +++ b/InnerDetector/InDetDetDescr/InDetGeoModelUtils/src/VolumeSplitterUtils.cxx @@ -1,6 +1,6 @@ /* - Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration -*/ + Copyright (C) 2002-2017 CERN for the benefit of the ATLAS collaboration + */ #include "InDetGeoModelUtils/VolumeSplitterUtils.h" @@ -9,449 +9,424 @@ #include <iostream> namespace InDetDD { + const SegmentList& + SegmentSplitter::split(const Zone* zone, const Ray& ray) { + addChildSegment(zone, ray); + return m_segments; + } -const SegmentList & -SegmentSplitter::split(const Zone * zone, const Ray & ray) -{ - addChildSegment(zone, ray); - return m_segments; -} - -Ray -SegmentSplitter::addChildSegment(const Zone * zone, const Ray & ray) -{ - //std::cout << "addChildSegment: " << zone->label() << " " << ray << std::endl; - // Point startPoint = ray.start(); - // If start point is not yet know to be in current zone. We need to search children. - // If not in children searchPoint will return a copy of ray. Otherwise it will - // add segments of the child and return a ray that starts from the exit of the child. - Ray nextRay = ray; - if (!ray.foundStart()) { - nextRay = searchPoint(zone, ray); - } - bool exit = false; - while (nextRay.valid() && !exit) { - Point nextPoint = getNextBoundary(zone, nextRay); - addSegment(zone, nextRay.start(), nextPoint); - - //std::cout << "Next boundary: " << nextPoint; - //if (nextPoint.last()) std::cout << " Last "; - //if (nextPoint.exit()) std::cout << " Exit "; - //if (nextPoint.child()) std::cout << " Child: " << nextPoint.child()->label(); - //std::cout << std::endl; - - if (nextPoint.last()) { // last indicates end of ray. ie ray ended before exit or child entry. - nextRay.setInvalid(); - } else { - nextRay.set(nextPoint, nextRay.end()); + Ray + SegmentSplitter::addChildSegment(const Zone* zone, const Ray& ray) { + //std::cout << "addChildSegment: " << zone->label() << " " << ray << std::endl; + // Point startPoint = ray.start(); + // If start point is not yet know to be in current zone. We need to search children. + // If not in children searchPoint will return a copy of ray. Otherwise it will + // add segments of the child and return a ray that starts from the exit of the child. + Ray nextRay = ray; + + if (!ray.foundStart()) { + nextRay = searchPoint(zone, ray); } - exit = nextPoint.exit(); - if (nextPoint.child()) { - nextRay = addChildSegment(nextPoint.child(), nextRay); + bool exit = false; + while (nextRay.valid() && !exit) { + Point nextPoint = getNextBoundary(zone, nextRay); + addSegment(zone, nextRay.start(), nextPoint); + + //std::cout << "Next boundary: " << nextPoint; + //if (nextPoint.last()) std::cout << " Last "; + //if (nextPoint.exit()) std::cout << " Exit "; + //if (nextPoint.child()) std::cout << " Child: " << nextPoint.child()->label(); + //std::cout << std::endl; + + if (nextPoint.last()) { // last indicates end of ray. ie ray ended before exit or child entry. + nextRay.setInvalid(); + } else { + nextRay.set(nextPoint, nextRay.end()); + } + exit = nextPoint.exit(); + if (nextPoint.child()) { + nextRay = addChildSegment(nextPoint.child(), nextRay); + } } + return nextRay; } - return nextRay; -} -void -SegmentSplitter::addSegment(const Zone * zone, const Point & start, const Point & end) -{ - m_segments.add(zone->label(), start, end, zone->rotated()); -} + void + SegmentSplitter::addSegment(const Zone* zone, const Point& start, const Point& end) { + m_segments.add(zone->label(), start, end, zone->rotated()); + } -Point -SegmentSplitter::getNextBoundary(const Zone * zone, const Ray & ray) { - Point nextPoint; //invalid - // Search children - for(Zone::ChildIterator iter = zone->begin(); iter!= zone->end(); ++iter) { - const Zone * child = *iter; - // return invalid if not crossed. - Point newPoint = child->findEntry(ray); - nextPoint = nearestPoint(newPoint, nextPoint); - } - //search exit - nextPoint = nearestPoint(zone->findExit(ray), nextPoint); - if (!nextPoint.valid()) { - // Return endpoint - return ray.end(); - } - return nextPoint; -} + Point + SegmentSplitter::getNextBoundary(const Zone* zone, const Ray& ray) { + Point nextPoint; //invalid -Point SegmentSplitter::nearestPoint(const Point & point1, const Point & point2) -{ - if (!point2.valid()) return point1; - if (!point1.valid()) return point2; - if ((point2.r() == point1.r() && point2.z() < point1.z()) || - (point2.z() == point1.z() && point2.r() < point1.r())) { - return point2; - } - return point1; - //return first point. If either invalid return the other. - // If both invalid return invalid -} + // Search children + for (Zone::ChildIterator iter = zone->begin(); iter != zone->end(); ++iter) { + const Zone* child = *iter; + // return invalid if not crossed. + Point newPoint = child->findEntry(ray); + nextPoint = nearestPoint(newPoint, nextPoint); + } + //search exit + nextPoint = nearestPoint(zone->findExit(ray), nextPoint); + if (!nextPoint.valid()) { + // Return endpoint + return ray.end(); + } + return nextPoint; + } + + Point + SegmentSplitter::nearestPoint(const Point& point1, const Point& point2) { + if (!point2.valid()) return point1; + if (!point1.valid()) return point2; -Ray -SegmentSplitter::searchPoint(const Zone * zone, const Ray & ray) -{ - // If not found in children return original ray. - //std::cout << "Searching for point " << ray.start() << std::endl; - Ray nextRay = ray; - for(Zone::ChildIterator iter = zone->begin(); iter != zone->end(); ++iter) { - const Zone * child = *iter; - if (child->inSide(ray.start())) { - //std::cout << "Point " << ray.start() << " found in zone " << child->label() << std::endl; - nextRay = addChildSegment(child, ray); - break; + if ((point2.r() == point1.r() && point2.z() < point1.z()) || + (point2.z() == point1.z() && point2.r() < point1.r())) { + return point2; } + return point1; + //return first point. If either invalid return the other. + // If both invalid return invalid } - nextRay.setFound(); - return nextRay; -} -Zone::Zone(const std::string & label, bool rotated) - : m_label(label), + Ray + SegmentSplitter::searchPoint(const Zone* zone, const Ray& ray) { + // If not found in children return original ray. + //std::cout << "Searching for point " << ray.start() << std::endl; + Ray nextRay = ray; + + for (Zone::ChildIterator iter = zone->begin(); iter != zone->end(); ++iter) { + const Zone* child = *iter; + if (child->inSide(ray.start())) { + //std::cout << "Point " << ray.start() << " found in zone " << child->label() << std::endl; + nextRay = addChildSegment(child, ray); + break; + } + } + nextRay.setFound(); + return nextRay; + } + + Zone::Zone(const std::string& label, bool rotated) + : m_label(label), m_rotated(rotated) -{} + {} -Zone::~Zone() -{ - for (ChildIterator iter = begin(); iter!=end(); ++iter) { - delete *iter; + Zone::~Zone() { + for (ChildIterator iter = begin(); iter != end(); ++iter) { + delete *iter; + } } -} -void -Zone::add(const Zone * zone) -{ - m_children.push_back(zone); -} + void + Zone::add(const Zone* zone) { + m_children.push_back(zone); + } -UnboundedZone::UnboundedZone(const std::string & label) - : Zone(label) -{} + UnboundedZone::UnboundedZone(const std::string& label) + : Zone(label) + {} -bool -UnboundedZone::inSide(const Point &) const -{ - return true; -} + bool + UnboundedZone::inSide(const Point&) const { + return true; + } -Point -UnboundedZone::findEntry(const Ray &) const -{ - // Will never be called. - return Point(); -} + Point + UnboundedZone::findEntry(const Ray&) const { + // Will never be called. + return Point(); + } -Point -UnboundedZone::findExit(const Ray &) const -{ - // Invalid means no exit point. - return Point(); -} + Point + UnboundedZone::findExit(const Ray&) const { + // Invalid means no exit point. + return Point(); + } -TubeZone::TubeZone(const std::string & label, double zmin, double zmax, double rmin, double rmax, bool rotated) - : Zone(label, rotated), + TubeZone::TubeZone(const std::string& label, double zmin, double zmax, double rmin, double rmax, bool rotated) + : Zone(label, rotated), m_zmin(zmin), m_zmax(zmax), m_rmin(rmin), m_rmax(rmax) -{} + {} -TubeZone::TubeZone(const std::string & label, const GeoTube * shape, double zOffset, bool rotated) - : Zone(label, rotated), + TubeZone::TubeZone(const std::string& label, const GeoTube* shape, double zOffset, bool rotated) + : Zone(label, rotated), m_zmin(zOffset - shape->getZHalfLength()), m_zmax(zOffset + shape->getZHalfLength()), m_rmin(shape->getRMin()), m_rmax(shape->getRMax()) -{} + {} - -bool -TubeZone::inSide(const Point & point) const -{ - return (inZ(point.z()) && inR(point.r())); -} -bool -TubeZone::inZ(double z) const { - return (z >= m_zmin && z < m_zmax); -} + bool + TubeZone::inSide(const Point& point) const { + return(inZ(point.z()) && inR(point.r())); + } -bool -TubeZone::inR(double r) const { - return (r >= m_rmin && r < m_rmax); -} + bool + TubeZone::inZ(double z) const { + return(z >= m_zmin && z < m_zmax); + } + + bool + TubeZone::inR(double r) const { + return(r >= m_rmin && r < m_rmax); + } // Assume either vertical or horizontal. -Point -TubeZone::findEntry(const Ray & ray) const -{ - if (ray.horizontal()) { - if (inR(ray.start().r()) && ray.start().z() < m_zmin && ray.end().z() > m_zmin) { - Point p(m_zmin, ray.start().r()); - p.setChild(this); - return p; + Point + TubeZone::findEntry(const Ray& ray) const { + if (ray.horizontal()) { + if (inR(ray.start().r()) && ray.start().z() < m_zmin && ray.end().z() > m_zmin) { + Point p(m_zmin, ray.start().r()); + p.setChild(this); + return p; + } + } else if (ray.vertical()) { + if (inZ(ray.start().z()) && ray.start().r() < m_rmin && ray.end().r() > m_rmin) { + Point p(ray.start().z(), m_rmin); + p.setChild(this); + return p; + } + } else { + std::cout << "Unexpected case" << std::endl; } - } else if(ray.vertical()) { - if (inZ(ray.start().z()) && ray.start().r() < m_rmin && ray.end().r() > m_rmin) { - Point p(ray.start().z(), m_rmin); - p.setChild(this); - return p; - } - } else { - std::cout << "Unexpected case" << std::endl; - } - // Return invalid point since doesn't intersect. - return Point(); // invalid point -} - - -// Assume already inside. -Point -TubeZone::findExit(const Ray & ray) const -{ - if (ray.horizontal()) { - if (ray.end().z() > m_zmax) { - Point p(m_zmax, ray.start().r()); - p.setExit(); - return p; + // Return invalid point since doesn't intersect. + return Point(); // invalid point + } + +// Assume already inside. + Point + TubeZone::findExit(const Ray& ray) const { + if (ray.horizontal()) { + if (ray.end().z() > m_zmax) { + Point p(m_zmax, ray.start().r()); + p.setExit(); + return p; + } + } else if (ray.vertical()) { + if (ray.end().r() > m_rmax) { + Point p(ray.start().z(), m_rmax); + p.setExit(); + return p; + } + } else { + std::cout << "Unexpected case" << std::endl; } - } else if(ray.vertical()) { - if (ray.end().r() > m_rmax) { - Point p(ray.start().z(), m_rmax); - p.setExit(); - return p; + // ends with. Return invalid point. + return Point(); + } + + PconZone::PconZone(const std::string& label, bool rotated) + : Zone(label, rotated) + {} + + PconZone::PconZone(const std::string& label, const GeoPcon* shape, bool rotated) + : Zone(label, rotated) { + for (unsigned int i = 0; i < shape->getNPlanes(); ++i) { + addPlane(shape->getZPlane(i), shape->getRMinPlane(i), shape->getRMaxPlane(i)); } - } else { - std::cout << "Unexpected case" << std::endl; } - // ends with. Return invalid point. - return Point(); -} -PconZone::PconZone(const std::string & label, bool rotated) - : Zone(label, rotated) -{} + void + PconZone::addPlane(double z, double rmin, double rmax) { + m_z.push_back(z); + m_rmin.push_back(rmin); + m_rmax.push_back(rmax); + } -PconZone::PconZone(const std::string & label, const GeoPcon * shape, bool rotated) - : Zone(label, rotated) -{ - for (unsigned int i = 0; i < shape->getNPlanes(); ++i) { - addPlane(shape->getZPlane(i), shape->getRMinPlane(i), shape->getRMaxPlane(i)); + bool + PconZone::inSide(const Point& point) const { + // Assume comes in pairs with same CLHEP::radii. + for (unsigned int i = 0; i + 1 < m_z.size(); i += 2) { + if (inZ(i, point.z()) && inR(i, point.r())) return true; + } + return false; } -} - -void -PconZone::addPlane(double z, double rmin, double rmax) -{ - m_z.push_back(z); - m_rmin.push_back(rmin); - m_rmax.push_back(rmax); -} -bool -PconZone::inSide(const Point & point) const -{ - // Assume comes in pairs with same CLHEP::radii. - for (unsigned int i = 0; i+1 < m_z.size(); i += 2) { - if (inZ(i,point.z()) && inR(i,point.r())) return true; + bool + PconZone::inZ(unsigned int i, double z) const { + return(z >= m_z[i] && z < m_z[i + 1]); } - return false; -} -bool -PconZone::inZ(unsigned int i,double z) const { - return (z >= m_z[i] && z < m_z[i+1]); -} + bool + PconZone::inR(unsigned int i, double r) const { + if (i >= m_z.size()) return false; -bool -PconZone::inR(unsigned int i, double r) const { - if (i >= m_z.size()) return false; - return (r >= m_rmin[i] && r < m_rmax[i]); -} + return(r >= m_rmin[i] && r < m_rmax[i]); + } // Assume either vertical or horizontal. -Point -PconZone::findEntry(const Ray & ray) const -{ - if (ray.horizontal()) { - for (unsigned int i = 0; i+1 < m_z.size(); i += 2) { - if (inR(i, ray.start().r()) && ray.start().z() < m_z[i] && ray.end().z() > m_z[i]) { - Point p(m_z[i], ray.start().r()); - p.setChild(this); - return p; + Point + PconZone::findEntry(const Ray& ray) const { + if (ray.horizontal()) { + for (unsigned int i = 0; i + 1 < m_z.size(); i += 2) { + if (inR(i, ray.start().r()) && ray.start().z() < m_z[i] && ray.end().z() > m_z[i]) { + Point p(m_z[i], ray.start().r()); + p.setChild(this); + return p; + } } - } - } else if(ray.vertical()) { - for (unsigned int i = 0; i+1 < m_z.size(); i += 2) { - if (inZ(i, ray.start().z()) && ray.start().r() < m_rmin[i] && ray.end().r() > m_rmin[i]) { - Point p(ray.start().z(), m_rmin[i]); - p.setChild(this); - return p; - } - } - } else { - std::cout << "Unexpected case" << std::endl; - } - // Return invalid point since doesn't intersect. - return Point(); // invalid point -} - - -// Assume already inside. -Point -PconZone::findExit(const Ray & ray) const -{ - if (ray.horizontal()) { - for (unsigned int i = 0; i+1 < m_z.size(); i += 2) { - if (inZ(i,ray.start().z()) && ray.end().z() > m_z[i+1] && !inR(i+2, ray.start().r())) { - Point p(m_z[i+1], ray.start().r()); - p.setExit(); - return p; + } else if (ray.vertical()) { + for (unsigned int i = 0; i + 1 < m_z.size(); i += 2) { + if (inZ(i, ray.start().z()) && ray.start().r() < m_rmin[i] && ray.end().r() > m_rmin[i]) { + Point p(ray.start().z(), m_rmin[i]); + p.setChild(this); + return p; + } } - } - } else if(ray.vertical()) { - for (unsigned int i = 0; i+1 < m_z.size(); i += 2) { - if (inZ(i,ray.start().z()) && ray.end().r() > m_rmax[i]) { - Point p(ray.start().z(), m_rmax[i]); - p.setExit(); - return p; + } else { + std::cout << "Unexpected case" << std::endl; + } + // Return invalid point since doesn't intersect. + return Point(); // invalid point + } + +// Assume already inside. + Point + PconZone::findExit(const Ray& ray) const { + if (ray.horizontal()) { + for (unsigned int i = 0; i + 1 < m_z.size(); i += 2) { + if (inZ(i, ray.start().z()) && ray.end().z() > m_z[i + 1] && !inR(i + 2, ray.start().r())) { + Point p(m_z[i + 1], ray.start().r()); + p.setExit(); + return p; + } + } + } else if (ray.vertical()) { + for (unsigned int i = 0; i + 1 < m_z.size(); i += 2) { + if (inZ(i, ray.start().z()) && ray.end().r() > m_rmax[i]) { + Point p(ray.start().z(), m_rmax[i]); + p.setExit(); + return p; + } } + } else { + std::cout << "Unexpected case" << std::endl; } - } else { - std::cout << "Unexpected case" << std::endl; + // ends with. Return invalid point. + return Point(); } - // ends with. Return invalid point. - return Point(); -} -Ray::Ray() - : m_valid(false), + Ray::Ray() + : m_valid(false), m_found(false), m_horizontal(false), m_vertical(false) -{} + {} -Ray::Ray(const Point & start, const Point & end) + Ray::Ray(const Point& start, const Point& end) : m_valid(true), - m_found(false), - m_start(start), - m_end(end) -{ - m_end.setLast(); - setDirection(); -} + m_found(false), + m_start(start), + m_end(end) { + m_end.setLast(); + setDirection(); + } -void -Ray::set(const Point & start, const Point & end) -{ - m_start = start; - m_end = end; - setDirection(); -} + void + Ray::set(const Point& start, const Point& end) { + m_start = start; + m_end = end; + setDirection(); + } -void -Ray::setDirection() -{ - m_vertical = (m_start.z() == m_end.z()); - m_horizontal = (m_start.r() == m_end.r()); - if (m_vertical && m_horizontal) { - m_vertical = false; - m_horizontal = false; + void + Ray::setDirection() { + m_vertical = (m_start.z() == m_end.z()); + m_horizontal = (m_start.r() == m_end.r()); + if (m_vertical && m_horizontal) { + m_vertical = false; + m_horizontal = false; + } } -} -Point::Point() - : m_valid(false), + Point::Point() + : m_valid(false), m_exit(false), m_last(false), m_z(0), m_r(0), m_child(0) -{} + {} -Point::Point(double z, double r) - : m_valid(true), + Point::Point(double z, double r) + : m_valid(true), m_exit(false), m_last(false), m_z(z), m_r(r), m_child(0) -{} + {} -Segment::Segment(const std::string & label, const Point & start, const Point & end, bool rotated) - : m_label(label), + Segment::Segment(const std::string& label, const Point& start, const Point& end, bool rotated) + : m_label(label), m_rotated(rotated), m_zmin(start.z()), m_zmax(end.z()), m_rmin(start.r()), m_rmax(end.r()) -{} - -void -SegmentList::add(const std::string & label, const Point & start, const Point & end, bool rotated) -{ - m_segments.push_back(Segment(label, start, end, rotated)); -} + {} -void -SegmentList::add(const Segment & segment) -{ - m_segments.push_back(segment); -} - + void + SegmentList::add(const std::string& label, const Point& start, const Point& end, bool rotated) { + m_segments.push_back(Segment(label, start, end, rotated)); + } -bool -SegmentList::horizontal() const -{ - if (size() > 0) { - const Segment & seg = m_segments[0]; - return (seg.rmin() == seg.rmax()); + void + SegmentList::add(const Segment& segment) { + m_segments.push_back(segment); } - // Not relevant if array empty. - return true; -} -void -SegmentList::print() const -{ - for (std::vector<Segment>::const_iterator iter = m_segments.begin(); iter != m_segments.end(); ++iter) { - iter->print(); + bool + SegmentList::horizontal() const { + if (size() > 0) { + const Segment& seg = m_segments[0]; + return(seg.rmin() == seg.rmax()); + } + // Not relevant if array empty. + return true; } -} -void -Segment::print() const -{ - std::cout << m_label << " " - << m_zmin << " " - << m_zmax << " " - << m_rmin << " " - << m_rmax - << std::endl; -} - -std::ostream & operator<<(std::ostream & os, const InDetDD::Ray & ray) { - if (!ray.valid()) { - os << "INVALID"; - } else { - os << ray.start() << " --> " << ray.end(); - } - return os; -} + void + SegmentList::print() const { + for (std::vector<Segment>::const_iterator iter = m_segments.begin(); iter != m_segments.end(); ++iter) { + iter->print(); + } + } -std::ostream & operator<<(std::ostream & os, const InDetDD::Point & point) { - if (!point.valid()) { - os << "INVALID"; - } else { - os << "(" << point.z() << ", " << point.r() << ")"; + void + Segment::print() const { + std::cout << m_label << " " + << m_zmin << " " + << m_zmax << " " + << m_rmin << " " + << m_rmax + << std::endl; } - return os; -} - -} + std::ostream& + operator << (std::ostream& os, const InDetDD::Ray& ray) { + if (!ray.valid()) { + os << "INVALID"; + } else { + os << ray.start() << " --> " << ray.end(); + } + return os; + } + std::ostream& + operator << (std::ostream& os, const InDetDD::Point& point) { + if (!point.valid()) { + os << "INVALID"; + } else { + os << "(" << point.z() << ", " << point.r() << ")"; + } + return os; + } +}