From 82ee3ab7979d0a122487ffa8c889d37291b62331 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Wed, 23 Dec 2015 19:33:45 +0100 Subject: [PATCH 1/6] Provide hexgonal geometry for HGCal --- CondFormats/GeometryObjects/BuildFile.xml | 1 + .../interface/HGCalParameters.h | 61 ++ .../src/T_EventSetup_HGCalParameters.cc | 4 + CondFormats/GeometryObjects/src/classes.h | 1 + .../GeometryObjects/src/classes_def.xml | 37 + .../python/GeometryExtended2023Dev_cff.py | 2 + .../ForwardDetId/interface/HGCalDetId.h | 6 +- DataFormats/ForwardDetId/src/HGCalDetId.cc | 20 +- .../CaloTopology/interface/HGCalTopology.h | 7 + Geometry/CaloTopology/src/HGCalTopology.cc | 79 +- .../test/testHGCalTopology_cfg.py | 1 + Geometry/HGCalCommonData/data/v1/hgcalEE.xml | 1 + .../HGCalCommonData/data/v1/hgcalHEgem.xml | 1 + .../HGCalCommonData/data/v1/hgcalHEsci.xml | 1 + .../HGCalCommonData/data/v1/hgcalHEsil.xml | 1 + Geometry/HGCalCommonData/data/v2/hgcalEE.xml | 1 + .../HGCalCommonData/data/v2/hgcalHEgem.xml | 1 + .../HGCalCommonData/data/v2/hgcalHEsci.xml | 1 + .../HGCalCommonData/data/v2/hgcalHEsil.xml | 1 + .../HGCalCommonData/data/v3/hgcalHEsci.xml | 1 + Geometry/HGCalCommonData/data/v4/hgcalEE.xml | 1 + .../HGCalCommonData/data/v4/hgcalHEsci.xml | 1 + .../HGCalCommonData/data/v4/hgcalHEsil.xml | 1 + Geometry/HGCalCommonData/data/v5/hgcalEE.xml | 1 + .../HGCalCommonData/data/v5/hgcalHEsci.xml | 1 + .../HGCalCommonData/data/v5/hgcalHEsil.xml | 1 + .../HGCalCommonData/data/v6/hgcalHEsil.xml | 1 + .../HGCalCommonData/data/v7/hgcalCons.xml | 611 ++++++++------ Geometry/HGCalCommonData/data/v7/hgcalEE.xml | 17 +- .../HGCalCommonData/data/v7/hgcalHEsil.xml | 17 +- .../HGCalCommonData/data/v7/hgcalwafer.xml | 4 + .../interface/HGCalDDDConstants.h | 93 +-- .../interface/HGCalGeomParameters.h | 64 ++ .../interface/HGCalGeometryMode.h | 29 + .../interface/HGCalParametersFromDD.h | 17 + .../HGCalCommonData/plugins/BuildFile.xml | 5 +- .../plugins/DDHGCalWaferAlgo.cc | 8 +- .../plugins/DDHGCalWaferAlgo.h | 1 + .../plugins/HGCalNumberingInitialization.cc | 28 +- .../plugins/HGCalParametersESModule.cc | 59 ++ .../hgcalParametersInitialization_cfi.py | 19 + .../hgcalV6ParametersInitialization_cfi.py | 14 + .../python/testGeometryExtended_cff.py | 2 + .../HGCalCommonData/python/testHGCXML_cfi.py | 3 +- .../HGCalCommonData/src/HGCalDDDConstants.cc | 743 ++++++++--------- .../src/HGCalGeomParameters.cc | 752 ++++++++++++++++++ .../HGCalCommonData/src/HGCalGeometryMode.cc | 8 + .../src/HGCalParametersFromDD.cc | 90 +++ .../HGCalCommonData/test/HGCGeometryTester.cc | 21 +- .../test/HGCalNumberingTester.cc | 139 +--- .../test/testHGCGeometry_cfg.py | 63 +- .../test/testHGCalNumbering_cfg.py | 27 +- .../interface/HGCalGeometryLoader.h | 11 + .../plugins/HGCalGeometryESProducer.cc | 2 +- Geometry/HGCalGeometry/src/HGCalGeometry.cc | 100 ++- .../HGCalGeometry/src/HGCalGeometryLoader.cc | 187 +++-- .../HGCalGeometry/test/HGCalGeometryTester.cc | 117 +-- .../testHGCalGeometryFromLocalDB_cfg.py | 17 +- .../test/python/testHGCalGeometry_cfg.py | 20 +- .../HGCalSimProducers/src/HGCDigitizer.cc | 117 +-- .../CaloTest/src/HGCalTestNumbering.cc | 28 +- SimG4CMS/Calo/interface/HGCNumberingScheme.h | 7 +- SimG4CMS/Calo/interface/HGCSD.h | 10 +- SimG4CMS/Calo/src/HGCNumberingScheme.cc | 79 +- SimG4CMS/Calo/src/HGCSD.cc | 52 +- 65 files changed, 2676 insertions(+), 1140 deletions(-) create mode 100644 CondFormats/GeometryObjects/interface/HGCalParameters.h create mode 100644 CondFormats/GeometryObjects/src/T_EventSetup_HGCalParameters.cc create mode 100644 Geometry/HGCalCommonData/interface/HGCalGeomParameters.h create mode 100644 Geometry/HGCalCommonData/interface/HGCalGeometryMode.h create mode 100644 Geometry/HGCalCommonData/interface/HGCalParametersFromDD.h create mode 100644 Geometry/HGCalCommonData/plugins/HGCalParametersESModule.cc create mode 100644 Geometry/HGCalCommonData/python/hgcalParametersInitialization_cfi.py create mode 100644 Geometry/HGCalCommonData/python/hgcalV6ParametersInitialization_cfi.py create mode 100644 Geometry/HGCalCommonData/src/HGCalGeomParameters.cc create mode 100644 Geometry/HGCalCommonData/src/HGCalGeometryMode.cc create mode 100644 Geometry/HGCalCommonData/src/HGCalParametersFromDD.cc diff --git a/CondFormats/GeometryObjects/BuildFile.xml b/CondFormats/GeometryObjects/BuildFile.xml index d5bdb550c994e..e12416f9cbd1c 100644 --- a/CondFormats/GeometryObjects/BuildFile.xml +++ b/CondFormats/GeometryObjects/BuildFile.xml @@ -5,6 +5,7 @@ + diff --git a/CondFormats/GeometryObjects/interface/HGCalParameters.h b/CondFormats/GeometryObjects/interface/HGCalParameters.h new file mode 100644 index 0000000000000..5650937e412e8 --- /dev/null +++ b/CondFormats/GeometryObjects/interface/HGCalParameters.h @@ -0,0 +1,61 @@ +#ifndef CondFormats_GeometryObjects_HGCalParameters_h +#define CondFormats_GeometryObjects_HGCalParameters_h + +#include "CondFormats/Serialization/interface/Serializable.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +#include +#include +#include +#include + +class HGCalParameters { + +public: + + HGCalParameters(const std::string& nam): name_(nam) { } + ~HGCalParameters( void ) { } + + struct hgtrap { + int lay; + float bl, tl, h, dz, alpha, cellSize; + }; + + struct hgtrform { + int zp, lay, sec, subsec; + bool used; + CLHEP::Hep3Vector h3v; + CLHEP::HepRotation hr; + }; + + std::string name_; + int nCells_; + int nSectors_; + std::vector cellSize_; + std::vector modules_; + std::vector moduler_; + std::vector trform_; + std::vector layer_; + std::vector layerIndex_; + std::vector layerGroup_; + std::vector cellFactor_; + std::vector depth_; + std::vector depthIndex_; + std::vector depthLayerF_; + std::vector zLayerHex_; + std::vector rMinLayHex_; + std::vector rMaxLayHex_; + std::vector waferCopy_; + std::vector waferTypeL_; + std::vector waferTypeT_; + std::vector waferPos_; + std::vector cellFine_; + std::vector cellCoarse_; + std::vector layerGroupM_; + std::vector layerGroupO_; + std::vector boundR_; + double waferR_; + int mode_; + +}; + +#endif diff --git a/CondFormats/GeometryObjects/src/T_EventSetup_HGCalParameters.cc b/CondFormats/GeometryObjects/src/T_EventSetup_HGCalParameters.cc new file mode 100644 index 0000000000000..568ecaecbef40 --- /dev/null +++ b/CondFormats/GeometryObjects/src/T_EventSetup_HGCalParameters.cc @@ -0,0 +1,4 @@ +#include "CondFormats/GeometryObjects/interface/HGCalParameters.h" +#include "FWCore/Utilities/interface/typelookup.h" + +TYPELOOKUP_DATA_REG(HGCalParameters); diff --git a/CondFormats/GeometryObjects/src/classes.h b/CondFormats/GeometryObjects/src/classes.h index 98bd0610b6410..afcc0c9d18035 100644 --- a/CondFormats/GeometryObjects/src/classes.h +++ b/CondFormats/GeometryObjects/src/classes.h @@ -6,4 +6,5 @@ #include "CondFormats/GeometryObjects/interface/PTrackerParameters.h" #include "CondFormats/GeometryObjects/interface/HcalParameters.h" #include "CondFormats/GeometryObjects/interface/PHcalParameters.h" +#include "CondFormats/GeometryObjects/interface/HGCalParameters.h" diff --git a/CondFormats/GeometryObjects/src/classes_def.xml b/CondFormats/GeometryObjects/src/classes_def.xml index d9fe783699247..aaed2bca3a0cf 100644 --- a/CondFormats/GeometryObjects/src/classes_def.xml +++ b/CondFormats/GeometryObjects/src/classes_def.xml @@ -89,6 +89,43 @@ + + diff --git a/Configuration/Geometry/python/GeometryExtended2023Dev_cff.py b/Configuration/Geometry/python/GeometryExtended2023Dev_cff.py index e191ed7a7c739..802fdf4fb29d8 100644 --- a/Configuration/Geometry/python/GeometryExtended2023Dev_cff.py +++ b/Configuration/Geometry/python/GeometryExtended2023Dev_cff.py @@ -4,3 +4,5 @@ from Geometry.TrackerNumberingBuilder.trackerNumberingGeometry_cfi import * from Geometry.HcalCommonData.hcalParameters_cfi import * from Geometry.HcalCommonData.hcalDDDSimConstants_cfi import * +from Geometry.HGCalCommonData.hgcalParametersInitialization_cfi import * +from Geometry.HGCalCommonData.hgcalNumberingInitialization_cfi import * diff --git a/DataFormats/ForwardDetId/interface/HGCalDetId.h b/DataFormats/ForwardDetId/interface/HGCalDetId.h index daabc9eba0a55..ce4416f0fff48 100644 --- a/DataFormats/ForwardDetId/interface/HGCalDetId.h +++ b/DataFormats/ForwardDetId/interface/HGCalDetId.h @@ -19,6 +19,7 @@ class HGCalDetId : public DetId { static const int kHGCalLayerMask = 0x1F; static const int kHGCalZsideOffset = 24; static const int kHGCalZsideMask = 0x1; + static const int kHGCalMaskCell = 0xFFFFFF00; /** Create a null cellid*/ HGCalDetId(); @@ -30,6 +31,9 @@ class HGCalDetId : public DetId { HGCalDetId(const DetId& id); /** Assignment from a generic cell id */ HGCalDetId& operator=(const DetId& id); + + /** Converter for a geometry cell id */ + HGCalDetId geometryCell () const {return id_&kHGCalMaskCell;} /// get the absolute value of the cell #'s in x and y int cell() const { return id_&kHGCalCellMask; } @@ -38,7 +42,7 @@ class HGCalDetId : public DetId { int wafer() const { return (id_>>kHGCalWaferOffset)&kHGCalWaferMask; } /// get the wafer type - int waferType() const { return (id_>>kHGCalWaferTypeOffset)&kHGCalWaferTypeMask; } + int waferType() const { return ((id_>>kHGCalWaferTypeOffset)&kHGCalWaferTypeMask ? 1 : -1); } /// get the layer # int layer() const { return (id_>>kHGCalLayerOffset)&kHGCalLayerMask; } diff --git a/DataFormats/ForwardDetId/src/HGCalDetId.cc b/DataFormats/ForwardDetId/src/HGCalDetId.cc index 7d23b2c264f2e..7587bb787afa6 100644 --- a/DataFormats/ForwardDetId/src/HGCalDetId.cc +++ b/DataFormats/ForwardDetId/src/HGCalDetId.cc @@ -10,19 +10,21 @@ HGCalDetId::HGCalDetId(uint32_t rawid) : DetId(rawid) { } HGCalDetId::HGCalDetId(ForwardSubdetector subdet, int zp, int lay, int wafertype, int wafer, int cell) : DetId(Forward,subdet) { - - if (cell>kHGCalCellMask || cell<0 || wafer>kHGCalWaferMask || wafer<0 || wafertype>kHGCalWaferTypeMask || wafertype<0 || lay>kHGCalLayerMask || lay<0) { + + if (wafertype < 0) wafertype = 0; + if (cell>kHGCalCellMask || cell<0 || wafer>kHGCalWaferMask || wafer<0 || wafertype>kHGCalWaferTypeMask || lay>kHGCalLayerMask || lay<0) { #ifdef DebugLog - std::cout << "[HGCalDetId] request for new id for layer=" << lay + std::cout << "[HGCalDetId] request for new id for" + << " layer=" << lay << ":" << kHGCalLayerMask << " @ zp=" << zp - << " wafer=" << wafer - << " waferType=" << wafertype - << " cell=" << cell + << " wafer=" << wafer << ":" << kHGCalWaferMask + << " waferType=" << wafertype << ":" << kHGCalWaferTypeMask + << " cell=" << cell << ":" << kHGCalCellMask << " for subdet=" << subdet << " has one or more fields out of bounds and will be reset" << std::endl; #endif - zp = lay = wafertype = wafer = cell=0; + zp = lay = wafertype = wafer = cell = 0; } id_ |= ((cell & kHGCalCellMask) << kHGCalCellOffset); id_ |= ((wafer & kHGCalWaferMask) << kHGCalWaferOffset); @@ -45,7 +47,7 @@ bool HGCalDetId::isValid(ForwardSubdetector subdet, int zp, int lay, int waferty bool ok = ((subdet == HGCEE || subdet == HGCHEF || subdet == HGCHEB) && (cell >= 0 && cell <= kHGCalCellMask) && (wafer >= 0 && wafer <= kHGCalWaferMask) && - (wafertype >= 0 || wafertype <= kHGCalWaferTypeMask) && + (wafertype <= kHGCalWaferTypeMask) && (lay >= 0 && lay <= kHGCalLayerMask) && (zp == -1 || zp == 1)); #ifdef DebugLog @@ -54,7 +56,7 @@ bool HGCalDetId::isValid(ForwardSubdetector subdet, int zp, int lay, int waferty << (subdet == HGCEE || subdet == HGCHEF || subdet == HGCHEB) << " Cell " << cell << ":" << (cell >= 0 && cell <= kHGCalCellMask) << " Wafer " << wafer << ":" << (wafer >= 0 && wafer <= kHGCalWaferMask) - << " WaferType " << wafertype << ":" << (wafertype >= 0 || wafertype <= kHGCalWaferTypeMask) + << " WaferType " << wafertype << ":" << (wafertype <= kHGCalWaferTypeMask) << " Layer " << lay << ":" << (lay >= 0 && lay <= kHGCalLayerMask) << " zp " << zp << ":" << (zp == -1 || zp == 1) << std::endl; #endif diff --git a/Geometry/CaloTopology/interface/HGCalTopology.h b/Geometry/CaloTopology/interface/HGCalTopology.h index 33544b636769e..92538027a4da5 100644 --- a/Geometry/CaloTopology/interface/HGCalTopology.h +++ b/Geometry/CaloTopology/interface/HGCalTopology.h @@ -4,8 +4,10 @@ #include "DataFormats/ForwardDetId/interface/ForwardSubdetector.h" #include "DataFormats/ForwardDetId/interface/HGCEEDetId.h" #include "DataFormats/ForwardDetId/interface/HGCHEDetId.h" +#include "DataFormats/ForwardDetId/interface/HGCalDetId.h" #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h" #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h" +#include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h" #include #include @@ -82,6 +84,9 @@ class HGCalTopology : public CaloSubdetectorTopology { return vNeighborsDetId; } + ///Geometry mode + int geomMode() const {return hdcons_.geomMode();} + ///Dense indexing virtual uint32_t detId2denseId(const DetId& id) const; virtual DetId denseId2detId(uint32_t denseId) const; @@ -93,6 +98,7 @@ class HGCalTopology : public CaloSubdetectorTopology { unsigned int totalModules() const {return kSizeForDenseIndexing;} unsigned int totalGeomModules() const {return (unsigned int)(2*kHGeomHalf_);} + unsigned int allGeomModules() const; const HGCalDDDConstants& dddConstants () const {return hdcons_;} @@ -101,6 +107,7 @@ class HGCalTopology : public CaloSubdetectorTopology { DetId offsetBy(const DetId startId, int nrStepsX, int nrStepsY) const; DetId switchZSide(const DetId startId) const; + /// Use subSector in square mode as wafer type in hexagon mode static const int subSectors_ = 2; struct DecodedDetId { diff --git a/Geometry/CaloTopology/src/HGCalTopology.cc b/Geometry/CaloTopology/src/HGCalTopology.cc index ba1ad69719ac2..0d155eac851dc 100644 --- a/Geometry/CaloTopology/src/HGCalTopology.cc +++ b/Geometry/CaloTopology/src/HGCalTopology.cc @@ -10,8 +10,13 @@ HGCalTopology::HGCalTopology(const HGCalDDDConstants& hdcons, sectors_ = hdcons_.sectors(); layers_ = hdcons_.layers(true); cells_ = hdcons_.maxCells(true); - kHGhalf_ = sectors_*layers_*subSectors_*cells_ ; - kHGeomHalf_ = (half_ ? (sectors_*layers_*subSectors_) : (sectors_*layers_)); + if (geomMode() == HGCalGeometryMode::Square) { + kHGhalf_ = sectors_*layers_*subSectors_*cells_ ; + kHGeomHalf_ = (half_ ? (sectors_*layers_*subSectors_) : (sectors_*layers_)); + } else { + kHGhalf_ = sectors_*layers_*cells_ ; + kHGeomHalf_ = sectors_*layers_; + } kSizeForDenseIndexing = (unsigned int)(2*kHGhalf_); #ifdef DebugLog std::cout << "HGCalTopology initialized for subdetector " << subdet_ @@ -22,13 +27,26 @@ HGCalTopology::HGCalTopology(const HGCalDDDConstants& hdcons, #endif } +unsigned int HGCalTopology::allGeomModules() const { + int n = (geomMode() == HGCalGeometryMode::Square) ? + (2*kHGeomHalf_) : (2*hdcons_.wafers()); + return (unsigned int)(n); +} + uint32_t HGCalTopology::detId2denseId(const DetId& id) const { HGCalTopology::DecodedDetId id_ = decode(id); int isubsec= (id_.iSubSec > 0) ? 1 : 0; - uint32_t idx = (uint32_t)((((id_.zside > 0) ? kHGhalf_ : 0) + - ((((id_.iCell-1)*layers_+id_.iLay-1)*sectors_+ - id_.iSec-1)*subSectors_+isubsec))); + uint32_t idx; + if (geomMode() == HGCalGeometryMode::Square) { + idx = (uint32_t)((((id_.zside > 0) ? kHGhalf_ : 0) + + ((((id_.iCell-1)*layers_+id_.iLay-1)*sectors_+ + id_.iSec-1)*subSectors_+isubsec))); + } else { + idx = (uint32_t)((((id_.zside > 0) ? kHGhalf_ : 0) + + ((((id_.iCell-1)*layers_+id_.iLay-1)*sectors_+ + id_.iSec)*subSectors_+isubsec))); + } return idx; } @@ -40,7 +58,11 @@ DetId HGCalTopology::denseId2detId(uint32_t hi) const { int di = ((int)(hi)%kHGhalf_); int iSubSec= (di%subSectors_); id_.iSubSec= (iSubSec == 0 ? -1 : 1); - id_.iSec = (((di-iSubSec)/subSectors_)%sectors_+1); + if (geomMode() == HGCalGeometryMode::Square) { + id_.iSec = (((di-iSubSec)/subSectors_)%sectors_+1); + } else { + id_.iSec = (((di-iSubSec)/subSectors_)%sectors_); + } id_.iLay = (((((di-iSubSec)/subSectors_)-id_.iSec+1)/sectors_)%layers_+1); id_.iCell = (((((di-iSubSec)/subSectors_)-id_.iSec+1)/sectors_-id_.iLay+1)/layers_+1); return encode(id_); @@ -53,18 +75,31 @@ uint32_t HGCalTopology::detId2denseGeomId(const DetId& id) const { HGCalTopology::DecodedDetId id_ = decode(id); int isubsec= (half_ && id_.iSubSec > 0) ? 1 : 0; - uint32_t idx = (uint32_t)(((id_.zside > 0) ? kHGeomHalf_ : 0) + - ((isubsec*layers_+id_.iLay-1)*sectors_+ - id_.iSec-1)); + uint32_t idx; + if (geomMode() == HGCalGeometryMode::Square) { + idx = (uint32_t)(((id_.zside > 0) ? kHGeomHalf_ : 0) + + ((isubsec*layers_+id_.iLay-1)*sectors_+id_.iSec-1)); + } else { + idx = (uint32_t)(((id_.zside > 0) ? kHGeomHalf_ : 0) + + ((isubsec*layers_+id_.iLay-1)*sectors_+id_.iSec)); + } return idx; } bool HGCalTopology::valid(const DetId& id) const { HGCalTopology::DecodedDetId id_ = decode(id); - bool flag = (id.det() == DetId::Forward && id.subdetId() == (int)(subdet_) && - id_.iCell >= 0 && id_.iCell < cells_ && id_.iLay > 0 && - id_.iLay <= layers_ && id_.iSec > 0 && id_.iSec <= sectors_); + bool flag; + if (geomMode() == HGCalGeometryMode::Square) { + flag = (id.det() == DetId::Forward && id.subdetId() == (int)(subdet_) && + id_.iCell >= 0 && id_.iCell < cells_ && id_.iLay > 0 && + id_.iLay <= layers_ && id_.iSec > 0 && id_.iSec <= sectors_); + } else { + flag = (id.det() == DetId::Forward && id.subdetId() == (int)(subdet_) && + id_.iCell >= 0 && id_.iCell < cells_ && id_.iLay > 0 && + id_.iLay <= layers_ && id_.iSec >= 0 && id_.iSec <= sectors_); + if (flag) flag = hdcons_.isValid(id_.iLay,id_.iSec,id_.iCell,true); + } return flag; } @@ -97,7 +132,11 @@ HGCalTopology::DecodedDetId HGCalTopology::geomDenseId2decId(const uint32_t& hi) int di = ((int)(hi)%kHGeomHalf_); int iSubSec= (di%subSectors_); id_.iSubSec= (iSubSec == 0 ? -1 : 1); - id_.iSec = (((di-iSubSec)/subSectors_)%sectors_+1); + if (geomMode() == HGCalGeometryMode::Square) { + id_.iSec = (((di-iSubSec)/subSectors_)%sectors_+1); + } else { + id_.iSec = (((di-iSubSec)/subSectors_)%sectors_); + } id_.iLay = (((((di-iSubSec)/subSectors_)-id_.iSec+1)/sectors_)%layers_+1); } return id_; @@ -106,7 +145,15 @@ HGCalTopology::DecodedDetId HGCalTopology::geomDenseId2decId(const uint32_t& hi) HGCalTopology::DecodedDetId HGCalTopology::decode(const DetId& startId) const { HGCalTopology::DecodedDetId id_; - if (subdet_ == HGCEE) { + if (geomMode() == HGCalGeometryMode::Hexagon) { + HGCalDetId id(startId); + id_.iCell = id.cell(); + id_.iLay = id.layer(); + id_.iSec = id.wafer(); + id_.iSubSec= id.waferType(); + id_.zside = id.zside(); + id_.subdet = id.subdetId(); + } else if (subdet_ == HGCEE) { HGCEEDetId id(startId); id_.iCell = id.cell(); id_.iLay = id.layer(); @@ -130,7 +177,9 @@ DetId HGCalTopology::encode(const HGCalTopology::DecodedDetId& id_) const { int isubsec= (id_.iSubSec > 0) ? 1 : 0; DetId id; - if (subdet_ == HGCEE) { + if (geomMode() == HGCalGeometryMode::Hexagon) { + id = HGCalDetId(subdet_,id_.zside,id_.iLay,isubsec,id_.iSec,id_.iCell).rawId(); + } else if (subdet_ == HGCEE) { id = HGCEEDetId(subdet_,id_.zside,id_.iLay,id_.iSec,isubsec,id_.iCell).rawId(); } else { id = HGCHEDetId(subdet_,id_.zside,id_.iLay,id_.iSec,isubsec,id_.iCell).rawId(); diff --git a/Geometry/CaloTopology/test/testHGCalTopology_cfg.py b/Geometry/CaloTopology/test/testHGCalTopology_cfg.py index 6eef3db3630c6..7fd6622a4ddbc 100644 --- a/Geometry/CaloTopology/test/testHGCalTopology_cfg.py +++ b/Geometry/CaloTopology/test/testHGCalTopology_cfg.py @@ -4,6 +4,7 @@ process.load("SimGeneral.HepPDTESSource.pdt_cfi") process.load("Geometry.HGCalCommonData.testHGCalXML_cfi") +process.load("Geometry.HGCalCommonData.hgcalParametersInitialization_cfi") process.load("Geometry.HGCalCommonData.hgcalNumberingInitialization_cfi") process.load("Geometry.CaloEventSetup.HGCalTopology_cfi") diff --git a/Geometry/HGCalCommonData/data/v1/hgcalEE.xml b/Geometry/HGCalCommonData/data/v1/hgcalEE.xml index 0e73fd3cbe354..e463c957e7f00 100644 --- a/Geometry/HGCalCommonData/data/v1/hgcalEE.xml +++ b/Geometry/HGCalCommonData/data/v1/hgcalEE.xml @@ -110,6 +110,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v1/hgcalHEgem.xml b/Geometry/HGCalCommonData/data/v1/hgcalHEgem.xml index 981880a5b2f7f..d0385de5fe299 100644 --- a/Geometry/HGCalCommonData/data/v1/hgcalHEgem.xml +++ b/Geometry/HGCalCommonData/data/v1/hgcalHEgem.xml @@ -81,6 +81,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v1/hgcalHEsci.xml b/Geometry/HGCalCommonData/data/v1/hgcalHEsci.xml index 535e84312c949..c3983250ecf10 100644 --- a/Geometry/HGCalCommonData/data/v1/hgcalHEsci.xml +++ b/Geometry/HGCalCommonData/data/v1/hgcalHEsci.xml @@ -78,6 +78,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v1/hgcalHEsil.xml b/Geometry/HGCalCommonData/data/v1/hgcalHEsil.xml index 544547cedbd9e..6b18064beaf80 100644 --- a/Geometry/HGCalCommonData/data/v1/hgcalHEsil.xml +++ b/Geometry/HGCalCommonData/data/v1/hgcalHEsil.xml @@ -86,6 +86,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v2/hgcalEE.xml b/Geometry/HGCalCommonData/data/v2/hgcalEE.xml index 198b55d8ffa6b..b9aa2f1bd62c1 100644 --- a/Geometry/HGCalCommonData/data/v2/hgcalEE.xml +++ b/Geometry/HGCalCommonData/data/v2/hgcalEE.xml @@ -113,6 +113,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v2/hgcalHEgem.xml b/Geometry/HGCalCommonData/data/v2/hgcalHEgem.xml index ab4afeab7dfd1..eba5aa610abd1 100644 --- a/Geometry/HGCalCommonData/data/v2/hgcalHEgem.xml +++ b/Geometry/HGCalCommonData/data/v2/hgcalHEgem.xml @@ -78,6 +78,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v2/hgcalHEsci.xml b/Geometry/HGCalCommonData/data/v2/hgcalHEsci.xml index d0af8a6c0465c..6ccd0d4a39f3d 100644 --- a/Geometry/HGCalCommonData/data/v2/hgcalHEsci.xml +++ b/Geometry/HGCalCommonData/data/v2/hgcalHEsci.xml @@ -78,6 +78,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v2/hgcalHEsil.xml b/Geometry/HGCalCommonData/data/v2/hgcalHEsil.xml index 9ea2f44f0ff8c..eef36f0aab677 100644 --- a/Geometry/HGCalCommonData/data/v2/hgcalHEsil.xml +++ b/Geometry/HGCalCommonData/data/v2/hgcalHEsil.xml @@ -93,6 +93,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v3/hgcalHEsci.xml b/Geometry/HGCalCommonData/data/v3/hgcalHEsci.xml index 0a022d3bda309..b848f9a94c563 100644 --- a/Geometry/HGCalCommonData/data/v3/hgcalHEsci.xml +++ b/Geometry/HGCalCommonData/data/v3/hgcalHEsci.xml @@ -78,6 +78,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v4/hgcalEE.xml b/Geometry/HGCalCommonData/data/v4/hgcalEE.xml index 3c9cafca08541..f9076fe15b8e4 100644 --- a/Geometry/HGCalCommonData/data/v4/hgcalEE.xml +++ b/Geometry/HGCalCommonData/data/v4/hgcalEE.xml @@ -113,6 +113,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v4/hgcalHEsci.xml b/Geometry/HGCalCommonData/data/v4/hgcalHEsci.xml index c6ef2884265a3..3cdf344caa06a 100644 --- a/Geometry/HGCalCommonData/data/v4/hgcalHEsci.xml +++ b/Geometry/HGCalCommonData/data/v4/hgcalHEsci.xml @@ -78,6 +78,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v4/hgcalHEsil.xml b/Geometry/HGCalCommonData/data/v4/hgcalHEsil.xml index bebae3bdc309c..5f3a6aa9f2071 100644 --- a/Geometry/HGCalCommonData/data/v4/hgcalHEsil.xml +++ b/Geometry/HGCalCommonData/data/v4/hgcalHEsil.xml @@ -87,6 +87,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v5/hgcalEE.xml b/Geometry/HGCalCommonData/data/v5/hgcalEE.xml index 4aaf7815ea4bb..6e0375d158be1 100644 --- a/Geometry/HGCalCommonData/data/v5/hgcalEE.xml +++ b/Geometry/HGCalCommonData/data/v5/hgcalEE.xml @@ -127,6 +127,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v5/hgcalHEsci.xml b/Geometry/HGCalCommonData/data/v5/hgcalHEsci.xml index 5934aae54fc3f..5ead9dd9c7050 100644 --- a/Geometry/HGCalCommonData/data/v5/hgcalHEsci.xml +++ b/Geometry/HGCalCommonData/data/v5/hgcalHEsci.xml @@ -80,6 +80,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v5/hgcalHEsil.xml b/Geometry/HGCalCommonData/data/v5/hgcalHEsil.xml index 2ce3e676e74af..4a7b527caa7c3 100644 --- a/Geometry/HGCalCommonData/data/v5/hgcalHEsil.xml +++ b/Geometry/HGCalCommonData/data/v5/hgcalHEsil.xml @@ -94,6 +94,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v6/hgcalHEsil.xml b/Geometry/HGCalCommonData/data/v6/hgcalHEsil.xml index 7e516d9cdb5fc..288f8c129d03b 100644 --- a/Geometry/HGCalCommonData/data/v6/hgcalHEsil.xml +++ b/Geometry/HGCalCommonData/data/v6/hgcalHEsil.xml @@ -94,6 +94,7 @@ + diff --git a/Geometry/HGCalCommonData/data/v7/hgcalCons.xml b/Geometry/HGCalCommonData/data/v7/hgcalCons.xml index 164885e15e451..2bb4ad9527ef5 100644 --- a/Geometry/HGCalCommonData/data/v7/hgcalCons.xml +++ b/Geometry/HGCalCommonData/data/v7/hgcalCons.xml @@ -2,253 +2,384 @@ + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/Geometry/HGCalCommonData/data/v7/hgcalEE.xml b/Geometry/HGCalCommonData/data/v7/hgcalEE.xml index 6f625d6efa4c9..747d185155976 100644 --- a/Geometry/HGCalCommonData/data/v7/hgcalEE.xml +++ b/Geometry/HGCalCommonData/data/v7/hgcalEE.xml @@ -2,7 +2,10 @@ - + + + + @@ -81,4 +84,16 @@ + + + + + + + + + + + + diff --git a/Geometry/HGCalCommonData/data/v7/hgcalHEsil.xml b/Geometry/HGCalCommonData/data/v7/hgcalHEsil.xml index 243dbc6f78227..d3d9f87580eb2 100644 --- a/Geometry/HGCalCommonData/data/v7/hgcalHEsil.xml +++ b/Geometry/HGCalCommonData/data/v7/hgcalHEsil.xml @@ -2,7 +2,10 @@ - + + + + @@ -65,4 +68,16 @@ + + + + + + + + + + + + diff --git a/Geometry/HGCalCommonData/data/v7/hgcalwafer.xml b/Geometry/HGCalCommonData/data/v7/hgcalwafer.xml index 2fb4d4884a233..0529e04dcc66d 100644 --- a/Geometry/HGCalCommonData/data/v7/hgcalwafer.xml +++ b/Geometry/HGCalCommonData/data/v7/hgcalwafer.xml @@ -8,6 +8,8 @@ + + @@ -70,6 +72,7 @@ + HGCalCellCoarse, HGCalCellCoarseHalf @@ -123,6 +126,7 @@ + HGCalCellFine, HGCalCellFineHalf diff --git a/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h b/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h index 22a0bccc4af3e..17f413022cd0a 100644 --- a/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h +++ b/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h @@ -14,84 +14,77 @@ #include #include #include -#include - +#include "CondFormats/GeometryObjects/interface/HGCalParameters.h" #include "DetectorDescription/Core/interface/DDsvalues.h" -class DDCompactView; -class DDFilteredView; - class HGCalDDDConstants { public: - struct hgtrap { - hgtrap(int lay0, float bl0, float tl0, float h0, float dz0, float alpha0): - lay(lay0),bl(bl0),tl(tl0),h(h0),dz(dz0),alpha(alpha0),cellSize(0) {} - int lay; - float bl, tl, h, dz, alpha, cellSize; - }; - struct hgtrform { - hgtrform(int zp0, int lay0, int sec0, int subsec0): zp(zp0), lay(lay0), sec(sec0), subsec(subsec0),used(false) {} - int zp, lay, sec, subsec; - bool used; - CLHEP::Hep3Vector h3v; - CLHEP::HepRotation hr; - }; - - - HGCalDDDConstants(const DDCompactView& cpv, std::string & name); + HGCalDDDConstants(const HGCalParameters* hp, const std::string name); ~HGCalDDDConstants(); std::pair assignCell(float x, float y, int lay, int subSec, bool reco) const; - std::pair assignCell(float x, float y, float h, float bl, float tl, - float alpha, float cellSize) const; + std::pair assignCellSquare(float x, float y, float h, float bl, + float tl, float alpha, + float cellSize) const; + std::pair assignCellHexagon(float x, float y) const; + double cellSizeHex(int type) const; std::pair findCell(int cell, int lay, int subSec, bool reco) const; - std::pair findCell(int cell, float h, float bl, float tl, - float alpha, float cellSize) const; + std::pair findCellSquare(int cell, float h, float bl, float tl, + float alpha, float cellSize) const; + int geomMode() const {return hgpar_->mode_;} bool isValid(int lay, int mod, int cell, bool reco) const; - unsigned int layers(bool reco) const {return (reco ? depthIndex.size() : layerIndex.size());} - std::pair locateCell(int cell, int lay, int subSec, + unsigned int layers(bool reco) const; + std::pair locateCell(int cell, int lay, int type, bool reco) const; + std::pair locateCellHex(int cell, int wafer, bool reco) const; int maxCells(bool reco) const; int maxCells(int lay, bool reco) const; - int maxCells(float h, float bl, float tl, float alpha, - float cellSize) const; + int maxCellsSquare(float h, float bl, float tl, float alpha, + float cellSize) const; int maxRows(int lay, bool reco) const; + int modules(int lay, bool reco) const; std::pair newCell(int cell, int layer, int sector, int subsector, int incrx, int incry, bool half) const; std::pair newCell(int cell, int layer, int subsector, int incrz, bool half) const; int newCell(int kx, int ky, int lay, int subSec) const; std::vector numberCells(int lay, bool reco) const; - std::vector numberCells(float h, float bl, float tl, float alpha, - float cellSize) const; - int sectors() const {return nSectors;} - std::pair simToReco(int cell, int layer, bool half) const; + std::vector numberCellsSquare(float h, float bl, float tl, + float alpha, float cellSize) const; + int numberCellsHexagon(int wafer) const; + int sectors() const {return hgpar_->nSectors_;} + std::pair simToReco(int cell, int layer, int mod, bool half) const; + int waferFromCopy(int copy) const; + bool waferInLayer(int wafer, int lay, bool reco) const; + GlobalPoint waferPosition(int wafer) const {return hgpar_->waferPos_.at(wafer); } + int wafers() const; + int waferToCopy(int wafer) const {return ((wafer>=0)&&(wafer< (int)(hgpar_->waferCopy_.size()))) ? hgpar_->waferCopy_[wafer] : (int)(hgpar_->waferCopy_.size());} + int waferTypeT(int wafer) const {return ((wafer>=0)&&(wafer<(int)(hgpar_->waferTypeT_.size()))) ? hgpar_->waferTypeT_[wafer] : 0;} + - std::vector::const_iterator getFirstModule(bool reco=false) const { return (reco ? moduler_.begin() : modules_.begin()); } - std::vector::const_iterator getLastModule(bool reco=false) const { return (reco ? moduler_.end() : modules_.end()); } - std::vector::const_iterator getFirstTrForm() const { return trform_.begin(); } - std::vector::const_iterator getLastTrForm() const { return trform_.end(); } + std::vector::const_iterator getFirstModule(bool reco=false) const {return (reco ? hgpar_->moduler_.begin() : hgpar_->modules_.begin());} + std::vector::const_iterator getLastModule(bool reco=false) const {return (reco ? hgpar_->moduler_.end() : hgpar_->modules_.end());} + std::vector::const_iterator getModule(int wafer) const; + std::vector::const_iterator getFirstTrForm() const {return hgpar_->trform_.begin();} + std::vector::const_iterator getLastTrForm() const {return hgpar_->trform_.end(); } - const std::vector & getModules( void ) const { return moduler_; } - const std::vector & getTrForms( void ) const { return trform_; } + const std::vector & getModules() const {return hgpar_->moduler_; } + const std::vector & getTrForms() const {return hgpar_->trform_; } private: - void initialize(const DDCompactView& cpv, std::string name); - void loadGeometry(const DDFilteredView& fv, const std::string& tag); - void loadSpecPars(const DDFilteredView& fv); - std::vector getDDDArray(const std::string &, - const DDsvalues_type &, int &) const; + int cellHex(double xx, double yy, const double& cellR, + const std::vector& pos) const; std::pair getIndex(int lay, bool reco) const; + void getParameterSquare(int lay, int subSec, bool reco, float& h, float& bl, + float& tl, float& alpha) const; + bool waferInLayer(int wafer, int lay) const; - int nCells, nSectors, nLayers; - std::vector cellSize_; - std::vector modules_, moduler_; - std::vector trform_; - std::vector layer_, layerIndex; - std::vector layerGroup_, cellFactor_, depth_, depthIndex; + const HGCalParameters* hgpar_; + const double tan30deg_; + double rmax_; }; #endif diff --git a/Geometry/HGCalCommonData/interface/HGCalGeomParameters.h b/Geometry/HGCalCommonData/interface/HGCalGeomParameters.h new file mode 100644 index 0000000000000..51177eb884ce9 --- /dev/null +++ b/Geometry/HGCalCommonData/interface/HGCalGeomParameters.h @@ -0,0 +1,64 @@ +#ifndef HGCalCommonData_HGCalGeomParameters_h +#define HGCalCommonData_HGCalGeomParameters_h + +/** \class HGCalGeomParameters + * + * this class extracts some geometry constants from CompactView + * to be used by Reco Geometry/Topology + * + * $Date: 2015/06/25 00:06:50 $ + * \author Sunanda Banerjee, Fermilab + * + */ + +#include +#include +#include + +#include "DetectorDescription/Core/interface/DDsvalues.h" +#include "DataFormats/GeometryVector/interface/GlobalPoint.h" +class DDCompactView; +class DDFilteredView; +class HGCalParameters; + +class HGCalGeomParameters { + +public: + + HGCalGeomParameters(); + ~HGCalGeomParameters(); + void loadGeometrySquare(const DDFilteredView&, HGCalParameters&, + const std::string&); + void loadGeometryHexagon(const DDFilteredView&, HGCalParameters&, + const std::string&, const DDCompactView*, + const std::string&, const std::string&); + void loadSpecParsSquare(const DDFilteredView&, HGCalParameters&); + void loadSpecParsHexagon(const DDFilteredView&, HGCalParameters&, + const DDCompactView*, const std::string&, + const std::string&); + +private: + std::vector getDDDArray(const std::string&, const DDsvalues_type&, + int&); + std::pair cellPosition(const std::map& wafers, + std::map::const_iterator& itrf, + unsigned int num, double rmax, + double ymax, double xx, double yy, + unsigned int ncells); + + struct layerParameters { + double rmin, rmax, zpos; + layerParameters(double rin=0, double rout=0, + double zp=0) : rmin(rin), rmax(rout), zpos(zp) {} + }; + struct cellParameters { + bool half; + GlobalPoint xyz; + cellParameters(bool h=false, + GlobalPoint p=GlobalPoint(0,0,0)) : half(h), xyz(p) {} + }; + + double waferSize_; +}; + +#endif diff --git a/Geometry/HGCalCommonData/interface/HGCalGeometryMode.h b/Geometry/HGCalCommonData/interface/HGCalGeometryMode.h new file mode 100644 index 0000000000000..61bcfe1cc6de0 --- /dev/null +++ b/Geometry/HGCalCommonData/interface/HGCalGeometryMode.h @@ -0,0 +1,29 @@ +#ifndef Geometry_HGCalCommonData_HGCalGeometryMode_H +#define Geometry_HGCalCommonData_HGCalGeometryMode_H + +#include "FWCore/Utilities/interface/Exception.h" +#include +#include +#include + +template< typename T > +class HGCalStringToEnumParser { + std::map enumMap; +public: + + HGCalStringToEnumParser(void); + + T parseString(const std::string &value) { + typename std::map::const_iterator itr = enumMap.find(value); + if (itr == enumMap.end()) + throw cms::Exception("Configuration") << "the value " << value + << " is not defined."; + return itr->second; + } +}; + +namespace HGCalGeometryMode { + enum GeometryMode { Square=0, Hexagon=1 }; +} + +#endif diff --git a/Geometry/HGCalCommonData/interface/HGCalParametersFromDD.h b/Geometry/HGCalCommonData/interface/HGCalParametersFromDD.h new file mode 100644 index 0000000000000..2651c5a5763f8 --- /dev/null +++ b/Geometry/HGCalCommonData/interface/HGCalParametersFromDD.h @@ -0,0 +1,17 @@ +#ifndef HGCalCommonData_HGCalParametersFromDD_h +#define HGCalCommonData_HGCalParametersFromDD_h + +#include +class DDCompactView; +class HGCalParameters; + +class HGCalParametersFromDD { +public: + HGCalParametersFromDD() {} + virtual ~HGCalParametersFromDD() {} + + bool build(const DDCompactView*, HGCalParameters&, const std::string&, + const std::string&, const std::string&); +}; + +#endif diff --git a/Geometry/HGCalCommonData/plugins/BuildFile.xml b/Geometry/HGCalCommonData/plugins/BuildFile.xml index 0be58bb4faf37..ef9852be73b3c 100644 --- a/Geometry/HGCalCommonData/plugins/BuildFile.xml +++ b/Geometry/HGCalCommonData/plugins/BuildFile.xml @@ -4,10 +4,11 @@ + - - + + diff --git a/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.cc b/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.cc index 7ed71d08334d9..e9f6a166e0897 100644 --- a/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.cc +++ b/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.cc @@ -13,7 +13,7 @@ #include "Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.h" #include "CLHEP/Units/GlobalSystemOfUnits.h" -//#define DebugLog +#define DebugLog DDHGCalWaferAlgo::DDHGCalWaferAlgo() { #ifdef DebugLog @@ -30,6 +30,7 @@ void DDHGCalWaferAlgo::initialize(const DDNumericArguments & nArgs, const DDStringVectorArguments & vsArgs) { cellSize = nArgs["CellSize"]; + cellType = (int)(nArgs["CellType"]); childNames = vsArgs["ChildNames"]; positionX = dbl_to_int(vArgs["PositionX"]); positionY = dbl_to_int(vArgs["PositionY"]); @@ -87,10 +88,11 @@ void DDHGCalWaferAlgo::execute(DDCompactView& cpv) { double xpos = dx*positionX[k]; double ypos = dy*positionY[k]; DDTranslation tran(xpos, ypos, 0); - cpv.position(DDName(name,idNameSpace), parentName, k, tran, rotation); + int copy = cellType*1000+k; + cpv.position(DDName(name,idNameSpace), parentName, copy, tran, rotation); #ifdef DebugLog edm::LogInfo("HGCalGeom") << "DDHGCalWaferAlgo: " - << DDName(name,idNameSpace) << " number " << k + << DDName(name,idNameSpace) << " number " << copy << " positioned in " << parentName << " at " << tran << " with " << rotation; #endif diff --git a/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.h b/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.h index 3e2fa0860fc22..cc360ed7d305c 100644 --- a/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.h +++ b/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.h @@ -24,6 +24,7 @@ class DDHGCalWaferAlgo : public DDAlgorithm { private: double cellSize; //Cell Size + int cellType; //Type (1 fine; 2 coarse) std::vector childNames; //Names of children std::vector positionX; //Position in X std::vector positionY; //Position in Y diff --git a/Geometry/HGCalCommonData/plugins/HGCalNumberingInitialization.cc b/Geometry/HGCalCommonData/plugins/HGCalNumberingInitialization.cc index 5c08aff74bdce..6e6238b31ca9e 100644 --- a/Geometry/HGCalCommonData/plugins/HGCalNumberingInitialization.cc +++ b/Geometry/HGCalCommonData/plugins/HGCalNumberingInitialization.cc @@ -13,8 +13,6 @@ // // Original Author: Sunanda Banerjee // Created: Tue Mar 21 16:40:29 PDT 2013 -// $Id: HGCalNumberingInitialization.cc,v 1.0 2013/12/24 12:47:41 sunanda Exp $ -// // @@ -25,11 +23,11 @@ // user include files #include #include -#include +#include +#include +#include -#include -#include -#include +#include #include #include @@ -51,29 +49,29 @@ class HGCalNumberingInitialization : public edm::ESProducer { }; HGCalNumberingInitialization::HGCalNumberingInitialization(const edm::ParameterSet& iConfig) : hgcalDDDConst_(0) { + name_ = iConfig.getUntrackedParameter("Name"); + edm::LogInfo("HGCalGeom") << "HGCalNumberingInitialization for " << name_; #ifdef DebugLog - std::cout <<"constructing HGCalNumberingInitialization for " << name_ << std::endl; + std::cout << "HGCalNumberingInitialization for " << name_ << std::endl; #endif setWhatProduced(this, name_); } - HGCalNumberingInitialization::~HGCalNumberingInitialization() {} // ------------ method called to produce the data ------------ HGCalNumberingInitialization::ReturnType HGCalNumberingInitialization::produce(const IdealGeometryRecord& iRecord) { -#ifdef DebugLog - std::cout << "in HGCalNumberingInitialization::produce" << std::endl; -#endif + + edm::LogInfo("HGCalGeom") << "in HGCalNumberingInitialization::produce"; if (hgcalDDDConst_ == 0) { - edm::ESTransientHandle pDD; - iRecord.get(pDD); - hgcalDDDConst_ = new HGCalDDDConstants(*pDD, name_); + edm::ESHandle pHGpar; + iRecord.get(name_, pHGpar); + hgcalDDDConst_ = new HGCalDDDConstants(&(*pHGpar), name_); } - return std::auto_ptr (hgcalDDDConst_) ; + return ReturnType(hgcalDDDConst_) ; } //define this as a plug-in diff --git a/Geometry/HGCalCommonData/plugins/HGCalParametersESModule.cc b/Geometry/HGCalCommonData/plugins/HGCalParametersESModule.cc new file mode 100644 index 0000000000000..f9d28539c1349 --- /dev/null +++ b/Geometry/HGCalCommonData/plugins/HGCalParametersESModule.cc @@ -0,0 +1,59 @@ +#include "FWCore/Framework/interface/ModuleFactory.h" +#include "FWCore/Framework/interface/ESProducer.h" +#include "FWCore/Framework/interface/ESTransientHandle.h" +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "CondFormats/GeometryObjects/interface/HGCalParameters.h" +#include "DetectorDescription/Core/interface/DDCompactView.h" +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "Geometry/HGCalCommonData/interface/HGCalParametersFromDD.h" + +#include + +//#define DeugLog + +class HGCalParametersESModule : public edm::ESProducer { +public: + HGCalParametersESModule( const edm::ParameterSet & ); + ~HGCalParametersESModule( void ); + + typedef boost::shared_ptr ReturnType; + + ReturnType produce( const IdealGeometryRecord&); + +private: + std::string name_, namew_, namec_; +}; + +HGCalParametersESModule::HGCalParametersESModule(const edm::ParameterSet& iC) { + + name_ = iC.getUntrackedParameter("Name"); + namew_ = iC.getUntrackedParameter("NameW"); + namec_ = iC.getUntrackedParameter("NameC"); + edm::LogInfo("HGCalGeom") << "HGCalParametersESModule for " << name_ << ":" + << namew_ << ":" << namec_; +#ifdef DebugLog + std::cout << "HGCalParametersESModule for " << name_ << ":" << namew_ << ":" + << namec_ << std::endl; +#endif + setWhatProduced(this, name_); +} + +HGCalParametersESModule::~HGCalParametersESModule() {} + +HGCalParametersESModule::ReturnType +HGCalParametersESModule::produce(const IdealGeometryRecord& iRecord) { + edm::LogInfo("HGCalGeom") + << "HGCalParametersESModule::produce(const IdealGeometryRecord& iRecord)"; + edm::ESTransientHandle cpv; + iRecord.get(cpv); + + HGCalParameters* ptp = new HGCalParameters(name_); + HGCalParametersFromDD builder; + builder.build(&(*cpv), *ptp, name_, namew_, namec_); + + return ReturnType(ptp) ; +} + +//define this as a plug-in +DEFINE_FWK_EVENTSETUP_MODULE(HGCalParametersESModule); diff --git a/Geometry/HGCalCommonData/python/hgcalParametersInitialization_cfi.py b/Geometry/HGCalCommonData/python/hgcalParametersInitialization_cfi.py new file mode 100644 index 0000000000000..97f031beda358 --- /dev/null +++ b/Geometry/HGCalCommonData/python/hgcalParametersInitialization_cfi.py @@ -0,0 +1,19 @@ +import FWCore.ParameterSet.Config as cms + +hgcalEEParametersInitialize = cms.ESProducer("HGCalParametersESModule", + Name = cms.untracked.string("HGCalEESensitive"), + NameW = cms.untracked.string("HGCalWafer"), + NameC = cms.untracked.string("HGCalCell") +) + +hgcalHESiParametersInitialize = cms.ESProducer("HGCalParametersESModule", + Name = cms.untracked.string("HGCalHESiliconSensitive"), + NameW = cms.untracked.string("HGCalWafer"), + NameC = cms.untracked.string("HGCalCell") +) + +hgcalHEScParametersInitialize = cms.ESProducer("HGCalParametersESModule", + Name = cms.untracked.string("HGCalHEScintillatorSensitive"), + NameW = cms.untracked.string("HGCalWafer"), + NameC = cms.untracked.string("HGCalCell") +) diff --git a/Geometry/HGCalCommonData/python/hgcalV6ParametersInitialization_cfi.py b/Geometry/HGCalCommonData/python/hgcalV6ParametersInitialization_cfi.py new file mode 100644 index 0000000000000..991435a016838 --- /dev/null +++ b/Geometry/HGCalCommonData/python/hgcalV6ParametersInitialization_cfi.py @@ -0,0 +1,14 @@ +import FWCore.ParameterSet.Config as cms + +hgcalEEParametersInitialize = cms.ESProducer("HGCalParametersESModule", + Name = cms.untracked.string("HGCalEESensitive"), + NameW = cms.untracked.string("HGCalWafer"), + NameC = cms.untracked.string("HGCalCell") +) + +hgcalHESiParametersInitialize = cms.ESProducer("HGCalParametersESModule", + Name = cms.untracked.string("HGCalHESiliconSensitive"), + NameW = cms.untracked.string("HGCalWafer"), + NameC = cms.untracked.string("HGCalCell") +) + diff --git a/Geometry/HGCalCommonData/python/testGeometryExtended_cff.py b/Geometry/HGCalCommonData/python/testGeometryExtended_cff.py index 694b3ae36322d..44d06e9be99e6 100644 --- a/Geometry/HGCalCommonData/python/testGeometryExtended_cff.py +++ b/Geometry/HGCalCommonData/python/testGeometryExtended_cff.py @@ -3,3 +3,5 @@ from Geometry.HGCalCommonData.testHGCalXML_cfi import * from Geometry.TrackerNumberingBuilder.trackerNumberingGeometry_cfi import * from Geometry.HcalCommonData.hcalDDConstants_cff import * +from Geometry.HGCalCommonData.hgcalParametersInitialization_cfi import * +from Geometry.HGCalCommonData.hgcalNumberingInitialization_cfi import * diff --git a/Geometry/HGCalCommonData/python/testHGCXML_cfi.py b/Geometry/HGCalCommonData/python/testHGCXML_cfi.py index 79e324131f968..e433f38026069 100644 --- a/Geometry/HGCalCommonData/python/testHGCXML_cfi.py +++ b/Geometry/HGCalCommonData/python/testHGCXML_cfi.py @@ -12,7 +12,8 @@ 'Geometry/HGCalCommonData/data/v7/hgcal.xml', 'Geometry/HGCalCommonData/data/v7/hgcalEE.xml', 'Geometry/HGCalCommonData/data/v7/hgcalHEsil.xml', - 'Geometry/HGCalCommonData/data/v7/hgcalwafer.xml'), + 'Geometry/HGCalCommonData/data/v7/hgcalwafer.xml', + 'Geometry/HGCalCommonData/data/v7/hgcalCons.xml'), rootNodeName = cms.string('cms:OCMS') ) diff --git a/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc b/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc index c7e52cef5ae34..b75e8a3884167 100644 --- a/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc +++ b/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc @@ -2,6 +2,7 @@ #include "FWCore/MessageLogger/interface/MessageLogger.h" #include "FWCore/Utilities/interface/Exception.h" +#include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h" #include "DetectorDescription/Base/interface/DDutils.h" #include "DetectorDescription/Core/interface/DDValue.h" @@ -13,53 +14,65 @@ //#define DebugLog -const double k_ScaleFromDDD = 0.1; const double k_horizontalShift = 1.0; +const double k_ScaleFromDDD = 0.1; + +HGCalDDDConstants::HGCalDDDConstants(const HGCalParameters* hp, + const std::string name) : hgpar_(hp), tan30deg_(std::tan(30.0*CLHEP::deg)) { -HGCalDDDConstants::HGCalDDDConstants(const DDCompactView& cpv, - std::string& nam) { - initialize(cpv, nam); + if (geomMode() == HGCalGeometryMode::Square) { + rmax_ = 0; + edm::LogInfo("HGCalGeom") << "HGCalDDDConstants initialized for " << name + << " with " << layers(false) << ":" <waferR_) * std::cos(30.0*CLHEP::deg); + edm::LogInfo("HGCalGeom") << "HGCalDDDConstants initialized for " << name + << " with " << layers(false) << ":" < HGCalDDDConstants::assignCell(float x, float y, int lay, int subSec, bool reco) const { - std::pair index = getIndex(lay, reco); + std::pair cellAssignment(-1,-1); int i = index.first; - if (i < 0) return std::pair(-1,-1); - float alpha, h, bl, tl; - - //set default values - std::pair cellAssignment( (x>0)?1:0, -1 ); - if (reco) { - h = moduler_[i].h; - bl = moduler_[i].bl; - tl = moduler_[i].tl; - alpha= moduler_[i].alpha; - if ((subSec>0 && alpha<0) || (subSec<=0 && alpha>0)) alpha = -alpha; - cellAssignment=assignCell(x, y, h, bl, tl, alpha, index.second); + if (i < 0) return cellAssignment; + if (geomMode() == HGCalGeometryMode::Square) { + float alpha, h, bl, tl; + getParameterSquare(i,subSec,reco,h,bl,tl,alpha); + cellAssignment = assignCellSquare(x, y, h, bl, tl, alpha, index.second); } else { - h = modules_[i].h; - bl = modules_[i].bl; - tl = modules_[i].tl; - alpha= modules_[i].alpha; - cellAssignment=assignCell(x, y, h, bl, tl, alpha, index.second); + float xx = (reco) ? x : k_ScaleFromDDD*x; + float yy = (reco) ? y : k_ScaleFromDDD*y; + cellAssignment = assignCellHexagon(xx,yy); } - return cellAssignment; } -std::pair HGCalDDDConstants::assignCell(float x, float y, float h, - float bl,float tl,float alpha, - float cellSize) const { +std::pair HGCalDDDConstants::assignCellSquare(float x, float y, + float h, float bl, + float tl, float alpha, + float cellSize) const { float a = (alpha==0) ? (2*h/(tl-bl)) : (h/(tl-bl)); float b = 2*h*bl/(tl-bl); @@ -68,24 +81,23 @@ std::pair HGCalDDDConstants::assignCell(float x, float y, float h, int phiSector = (x0 > 0) ? 1 : 0; if (alpha < 0) {x0 -= 0.5*(tl+bl); phiSector = 0;} else if (alpha > 0) {x0 += 0.5*(tl+bl); phiSector = 1;} - - - //determine the i-y + + //determine the i-y int ky = floor((y+h)/cellSize); - if( ky*cellSize> y+h) ky--; - if(ky<0) return std::pair(-1,-1); - if( (ky+1)*cellSize < (y+h) ) ky++; + if (ky*cellSize> y+h) ky--; + if (ky<0) return std::pair(-1,-1); + if ((ky+1)*cellSize < (y+h) ) ky++; int max_ky_allowed=floor(2*h/cellSize); - if(ky>max_ky_allowed-1) return std::pair(-1,-1); - + if (ky>max_ky_allowed-1) return std::pair(-1,-1); + //determine the i-x //notice we substitute y by the top of the candidate cell to reduce the dead zones int kx = floor(fabs(x0)/cellSize); - if( kx*cellSize > fabs(x0) ) kx--; - if(kx<0) return std::pair(-1,-1); - if( (kx+1)*cellSize < fabs(x0) ) kx++; + if (kx*cellSize > fabs(x0) ) kx--; + if (kx<0) return std::pair(-1,-1); + if ((kx+1)*cellSize < fabs(x0)) kx++; int max_kx_allowed=floor( ((ky+1)*cellSize+b+k_horizontalShift*cellSize)/(a*cellSize) ); - if(kx>max_kx_allowed-1) return std::pair(-1,-1); + if (kx>max_kx_allowed-1) return std::pair(-1,-1); //count cells summing in rows until required height //notice the bottom of the cell must be used @@ -100,35 +112,55 @@ std::pair HGCalDDDConstants::assignCell(float x, float y, float h, return std::pair(phiSector,icell); } +std::pair HGCalDDDConstants::assignCellHexagon(float x, + float y) const { + double xx(x), yy(y); + //First the wafer + int wafer = cellHex(xx, yy, rmax_, hgpar_->waferPos_); + // Now the cell + xx -= hgpar_->waferPos_[wafer].x(); + yy -= hgpar_->waferPos_[wafer].y(); + int cell(0); + if (hgpar_->waferTypeT_[wafer] == 1) + cell = cellHex(xx, yy, 0.5*k_ScaleFromDDD*hgpar_->cellSize_[0], hgpar_->cellFine_); + else + cell = cellHex(xx, yy, 0.5*k_ScaleFromDDD*hgpar_->cellSize_[1], hgpar_->cellCoarse_); + return std::pair(wafer,cell); +} + +double HGCalDDDConstants::cellSizeHex(int type) const { + int indx = (type == 1) ? 0 : 1; + double cell = (0.5*k_ScaleFromDDD*hgpar_->cellSize_[indx]); + return cell; +} + +unsigned int HGCalDDDConstants::layers(bool reco) const { + return (reco ? hgpar_->depthIndex_.size() : hgpar_->layerIndex_.size()); +} + std::pair HGCalDDDConstants::findCell(int cell, int lay, int subSec, bool reco) const { std::pair index = getIndex(lay, reco); int i = index.first; if (i < 0) return std::pair(-1,-1); - float alpha, h, bl, tl; - if (reco) { - h = moduler_[i].h; - bl = moduler_[i].bl; - tl = moduler_[i].tl; - alpha= moduler_[i].alpha; - if ((subSec>0 && alpha<0) || (subSec<=0 && alpha>0)) alpha = -alpha; + if (geomMode() == HGCalGeometryMode::Hexagon) { + return std::pair(-1,-1); } else { - h = modules_[i].h; - bl = modules_[i].bl; - tl = modules_[i].tl; - alpha= modules_[i].alpha; + float alpha, h, bl, tl; + getParameterSquare(i,subSec,reco,h,bl,tl,alpha); + return findCellSquare(cell, h, bl, tl, alpha, index.second); } - return findCell(cell, h, bl, tl, alpha, index.second); } -std::pair HGCalDDDConstants::findCell(int cell, float h, float bl, - float tl, float alpha, - float cellSize) const { +std::pair HGCalDDDConstants::findCellSquare(int cell, float h, + float bl, float tl, + float alpha, + float cellSize) const { //check if cell number is meaningful if(cell<0) return std::pair(-1,-1); - + //parameterization of the boundary of the trapezoid float a = (alpha==0) ? (2*h/(tl-bl)) : (h/(tl-bl)); float b = 2*h*bl/(tl-bl); @@ -139,7 +171,7 @@ std::pair HGCalDDDConstants::findCell(int cell, float h, float bl, //check if adding all the cells in this row is above the required cell //notice the top of the cell is used to maximize space - int cellsInRow( floor( ((iky+1)*cellSize+b+k_horizontalShift*cellSize)/(a*cellSize) ) ); + int cellsInRow(floor(((iky+1)*cellSize+b+k_horizontalShift*cellSize)/(a*cellSize))); if (testCell+cellsInRow > cell) break; testCell += cellsInRow; ky++; @@ -149,48 +181,100 @@ std::pair HGCalDDDConstants::findCell(int cell, float h, float bl, return std::pair(kx,ky); } -bool HGCalDDDConstants::isValid(int lay, int mod, int cell, bool reco) const { +std::vector::const_iterator +HGCalDDDConstants::getModule(int wafer) const { + + std::vector::const_iterator itr; + int type = ((wafer>=0)&&(wafer<(int)(hgpar_->waferTypeL_.size()))) ? + hgpar_->waferTypeL_[wafer] : 3; + int indx(1); + for (itr = hgpar_->moduler_.begin(); itr != hgpar_->moduler_.end(); + ++itr, ++indx) { + if (type == indx) return itr; + } + return hgpar_->moduler_.end(); +} - bool ok = ((lay > 0 && lay <= (int)(layers(reco))) && - (mod > 0 && mod <= sectors()) && - (cell >=0 && cell <= maxCells(lay,reco))); +bool HGCalDDDConstants::isValid(int lay, int mod, int cell, bool reco) const { + bool ok(false); + int cellmax(0), modmax(0); + if (geomMode() == HGCalGeometryMode::Square) { + cellmax = maxCells(lay,reco); + modmax = sectors(); + ok = ((lay > 0 && lay <= (int)(layers(reco))) && + (mod > 0 && mod <= modmax) && + (cell >=0 && cell <= cellmax)); + } else { + modmax = modules(lay,reco); + ok = ((lay > 0 && lay <= (int)(layers(reco))) && + (mod > 0 && mod <= modmax)); + if (ok) { + cellmax = (hgpar_->waferTypeT_[mod]==1) ? + (int)(hgpar_->cellFine_.size()) : (int)(hgpar_->cellCoarse_.size()); + ok = (cell >=0 && cell <= cellmax); + } + } + #ifdef DebugLog if (!ok) std::cout << "HGCalDDDConstants: Layer " << lay << ":" << (lay > 0 && (lay <= (int)(layers(reco)))) << " Module " - << mod << ":" << (mod > 0 && mod <= sectors()) << " Cell " - << cell << ":" << (cell >=0 && cell <= maxCells(lay,reco)) + << mod << ":" << (mod > 0 && mod <= modmax) << " Cell " + << cell << ":" << (cell >=0 && cell <= cellmax) << ":" << maxCells(reco) << std::endl; #endif return ok; } std::pair HGCalDDDConstants::locateCell(int cell, int lay, - int subSec, bool reco) const { - + int type, bool reco) const { + // type refers to subsector # for square cell and wafer # for hexagon cell + float x(999999.), y(999999.); std::pair index = getIndex(lay, reco); int i = index.first; - if (i < 0) return std::pair(999999.,999999.); - std::pair kxy = findCell(cell, lay, subSec, reco); - float alpha, h, bl, tl; - if (reco) { - h = moduler_[i].h; - bl = moduler_[i].bl; - tl = moduler_[i].tl; - alpha= moduler_[i].alpha; - if ((subSec>0 && alpha<0) || (subSec<=0 && alpha>0)) alpha = -alpha; + if (i < 0) return std::pair(x,y); + if (geomMode() == HGCalGeometryMode::Square) { + std::pair kxy = findCell(cell, lay, type, reco); + float alpha, h, bl, tl; + getParameterSquare(i,type,reco,h,bl,tl,alpha); + float cellSize = index.second; + x = (kxy.first+0.5)*cellSize; + if (alpha < 0) x -= 0.5*(tl+bl); + else if (alpha > 0) x -= 0.5*(tl+bl); + if (type != 1) x = -x; + y = ((kxy.second+0.5)*cellSize-h); } else { - h = modules_[i].h; - bl = modules_[i].bl; - tl = modules_[i].tl; - alpha= modules_[i].alpha; + x = hgpar_->waferPos_[type].x(); + y = hgpar_->waferPos_[type].y(); + if (hgpar_->waferTypeT_[type] == 1) { + x += hgpar_->cellFine_[cell].x(); + y += hgpar_->cellFine_[cell].y(); + } else { + x += hgpar_->cellCoarse_[cell].x(); + y += hgpar_->cellCoarse_[cell].y(); + } + if (!reco) { + x /= k_ScaleFromDDD; + y /= k_ScaleFromDDD; + } + } + return std::pair(x,y); +} + +std::pair HGCalDDDConstants::locateCellHex(int cell, int wafer, + bool reco) const { + float x(0), y(0); + if (hgpar_->waferTypeT_[wafer] == 1) { + x = hgpar_->cellFine_[cell].x(); + y = hgpar_->cellFine_[cell].y(); + } else { + x = hgpar_->cellCoarse_[cell].x(); + y = hgpar_->cellCoarse_[cell].y(); + } + if (!reco) { + x /= k_ScaleFromDDD; + y /= k_ScaleFromDDD; } - float cellSize = index.second; - float x = (kxy.first+0.5)*cellSize; - if (alpha < 0) x -= 0.5*(tl+bl); - else if (alpha > 0) x -= 0.5*(tl+bl); - if (subSec != 1) x = -x; - float y = ((kxy.second+0.5)*cellSize-h); return std::pair(x,y); } @@ -198,7 +282,7 @@ int HGCalDDDConstants::maxCells(bool reco) const { int cells(0); for (unsigned int i = 0; idepth_[i] : hgpar_->layer_[i]; if (cells < maxCells(lay, reco)) cells = maxCells(lay, reco); } return cells; @@ -209,23 +293,25 @@ int HGCalDDDConstants::maxCells(int lay, bool reco) const { std::pair index = getIndex(lay, reco); int i = index.first; if (i < 0) return 0; - float h, bl, tl, alpha; - if (reco) { - h = moduler_[i].h; - bl = moduler_[i].bl; - tl = moduler_[i].tl; - alpha= moduler_[i].alpha; + if (geomMode() == HGCalGeometryMode::Square) { + float h, bl, tl, alpha; + getParameterSquare(i,0,reco,h,bl,tl,alpha); + return maxCellsSquare(h, bl, tl, alpha, index.second); } else { - h = modules_[i].h; - bl = modules_[i].bl; - tl = modules_[i].tl; - alpha= modules_[i].alpha; + unsigned int cells(0); + for (unsigned int k=0; kwaferPos_.size(); ++k) { + if (waferInLayer(k,index.first)) { + unsigned int cell = (hgpar_->waferTypeT_[k]==1) ? + (hgpar_->cellFine_.size()) : (hgpar_->cellCoarse_.size()); + if (cell > cells) cells = cell; + } + } + return (int)(cells); } - return maxCells(h, bl, tl, alpha, index.second); } -int HGCalDDDConstants::maxCells(float h, float bl, float tl, float alpha, - float cellSize) const { +int HGCalDDDConstants::maxCellsSquare(float h, float bl, float tl, + float alpha, float cellSize) const { float a = (alpha==0) ? (2*h/(tl-bl)) : (h/(tl-bl)); float b = 2*h*bl/(tl-bl); @@ -233,25 +319,44 @@ int HGCalDDDConstants::maxCells(float h, float bl, float tl, float alpha, int ncells(0); //always use the top of the cell to maximize space int kymax = floor((2*h)/cellSize); - for (int iky=0; iky index = getIndex(lay, reco); int i = index.first; - if (i < 0) return 0; - float h = (reco) ? moduler_[i].h : modules_[i].h; - int kymax = floor((2*h)/index.second); + if (i < 0) return kymax; + if (geomMode() == HGCalGeometryMode::Square) { + float h = (reco) ? hgpar_->moduler_[i].h : hgpar_->modules_[i].h; + kymax = floor((2*h)/index.second); + } else { + for (unsigned int k=0; kwaferPos_.size(); ++k) { + if (waferInLayer(k,i)) { + int ky = ((hgpar_->waferCopy_[k])/100)%100; + if (ky > kymax) kymax = ky; + } + } + } return kymax; } +int HGCalDDDConstants::modules(int lay, bool reco) const { + int nmod(0); + std::pair index = getIndex(lay, reco); + if (index.first < 0) return nmod; + for (unsigned int k=0; kwaferPos_.size(); ++k) { + if (waferInLayer(k,index.first)) ++nmod; + } + return nmod; +} + std::pair HGCalDDDConstants::newCell(int cell, int layer, int sector, int subsector, int incrx, int incry, bool half) const { @@ -270,8 +375,8 @@ std::pair HGCalDDDConstants::newCell(int cell, int layer, int sector, kx -= maxCells(layer, true); sector += subsector; subsector =-subsector; - if (sector < 1) sector = nSectors; - else if (sector > nSectors) sector = 1; + if (sector < 1) sector = hgpar_->nSectors_; + else if (sector > hgpar_->nSectors_) sector = 1; } cell = newCell(kx, ky, layer, subSec); return std::pair(cell,sector*subsector); @@ -294,13 +399,13 @@ int HGCalDDDConstants::newCell(int kx, int ky, int lay, int subSec) const { std::pair index = getIndex(lay, true); int i = index.first; if (i < 0) return maxCells(true); - float alpha = (subSec == 0) ? modules_[i].alpha : subSec; + float alpha = (subSec == 0) ? hgpar_->modules_[i].alpha : subSec; float cellSize = index.second; float a = (alpha==0) ? - (2*moduler_[i].h/(moduler_[i].tl-moduler_[i].bl)) : - (moduler_[i].h/(moduler_[i].tl-moduler_[i].bl)); - float b = 2*moduler_[i].h*moduler_[i].bl/ - (moduler_[i].tl-moduler_[i].bl); + (2*hgpar_->moduler_[i].h/(hgpar_->moduler_[i].tl-hgpar_->moduler_[i].bl)) : + (hgpar_->moduler_[i].h/(hgpar_->moduler_[i].tl-hgpar_->moduler_[i].bl)); + float b = 2*hgpar_->moduler_[i].h*hgpar_->moduler_[i].bl/ + (hgpar_->moduler_[i].tl-hgpar_->moduler_[i].bl); int icell(kx); for (int iky=0; iky HGCalDDDConstants::numberCells(int lay, bool reco) const { std::pair index = getIndex(lay, reco); int i = index.first; + std::vector ncell; if (i >= 0) { - float h, bl, tl, alpha; - if (reco) { - h = moduler_[i].h; - bl = moduler_[i].bl; - tl = moduler_[i].tl; - alpha= moduler_[i].alpha; + if (geomMode() == HGCalGeometryMode::Square) { + float h, bl, tl, alpha; + getParameterSquare(i,0,reco,h,bl,tl,alpha); + return numberCellsSquare(h, bl, tl, alpha, index.second); } else { - h = modules_[i].h; - bl = modules_[i].bl; - tl = modules_[i].tl; - alpha= modules_[i].alpha; + for (unsigned int k=0; kwaferPos_.size(); ++k) { + if (waferInLayer(k,i)) { + unsigned int cell = (hgpar_->waferTypeT_[k]==1) ? + (hgpar_->cellFine_.size()) : (hgpar_->cellCoarse_.size()); + ncell.push_back((int)(cell)); + } + } } - return numberCells(h, bl, tl, alpha, index.second); - } else { - std::vector ncell; - return ncell; } + return ncell; } -std::vector HGCalDDDConstants::numberCells(float h, float bl, - float tl, float alpha, - float cellSize) const { +std::vector HGCalDDDConstants::numberCellsSquare(float h, float bl, + float tl, float alpha, + float cellSize) const { float a = (alpha==0) ? (2*h/(tl-bl)) : (h/(tl-bl)); float b = 2*h*bl/(tl-bl); @@ -344,299 +448,148 @@ std::vector HGCalDDDConstants::numberCells(float h, float bl, return ncell; } -std::pair HGCalDDDConstants::simToReco(int cell, int lay, +int HGCalDDDConstants::numberCellsHexagon(int wafer) const { + + int ncell(0); + if (wafer >= 0 && wafer < (int)(hgpar_->waferTypeT_.size())) { + if (hgpar_->waferTypeT_[wafer]==1) + ncell = (int)(hgpar_->cellFine_.size()); + else + ncell = (int)(hgpar_->cellCoarse_.size()); + } + return ncell; +} + +std::pair HGCalDDDConstants::simToReco(int cell, int lay, int mod, bool half) const { std::pair index = getIndex(lay, false); int i = index.first; if (i < 0) return std::pair(-1,-1); - float h = modules_[i].h; - float bl = modules_[i].bl; - float tl = modules_[i].tl; - float cellSize = cellFactor_[i]*index.second; - - std::pair kxy = findCell(cell, h, bl, tl, modules_[i].alpha, index.second); - int depth = layerGroup_[i]; - if(depth<0) return std::pair(-1,-1); - int kx = kxy.first/cellFactor_[i]; - int ky = kxy.second/cellFactor_[i]; - - float a = (half) ? (h/(tl-bl)) : (2*h/(tl-bl)); - float b = 2*h*bl/(tl-bl); - for (int iky=0; ikymodules_[i].h; + float bl = hgpar_->modules_[i].bl; + float tl = hgpar_->modules_[i].tl; + float cellSize = hgpar_->cellFactor_[i]*index.second; + std::pair kxy = findCellSquare(cell, h, bl, tl, hgpar_->modules_[i].alpha, index.second); + depth = hgpar_->layerGroup_[i]; + if (depth<0) return std::pair(-1,-1); + kx = kxy.first/hgpar_->cellFactor_[i]; + int ky = kxy.second/hgpar_->cellFactor_[i]; + + float a = (half) ? (h/(tl-bl)) : (2*h/(tl-bl)); + float b = 2*h*bl/(tl-bl); + for (int iky=0; ikywaferCopy_[k]) { + wafer = k; + break; + } } + return wafer; } -void HGCalDDDConstants::loadGeometry(const DDFilteredView& _fv, - const std::string & sdTag) { - - DDFilteredView fv = _fv; - bool dodet(true), first(true); - int zpFirst(0); - std::vector trforms; - - while (dodet) { - // DDTranslation t = fv.translation(); - const DDSolid & sol = fv.logicalPart().solid(); - std::string name = sol.name(); - int isd = (name.find(sdTag) == std::string::npos) ? -1 : 1; - if (isd > 0) { - std::vector copy = fv.copyNumbers(); - int nsiz = (int)(copy.size()); - int lay = (nsiz > 0) ? copy[nsiz-1] : -1; - int sec = (nsiz > 1) ? copy[nsiz-2] : -1; - int zp = (nsiz > 3) ? copy[nsiz-4] : -1; - if (zp !=1 ) zp = -1; - if (first) {first = false; zpFirst = zp;} - const DDTrap & trp = static_cast(sol); - HGCalDDDConstants::hgtrap mytr(lay,trp.x1(),trp.x2(), - 0.5*(trp.y1()+trp.y2()), - trp.halfZ(),trp.alpha1()); - int subs = (trp.alpha1()>0 ? 1 : 0); - if (std::find(layer_.begin(),layer_.end(),lay) == layer_.end()) { - for (unsigned int k=0; k indx = getIndex(lay, reco); + if (indx.first < 0) return false; + return waferInLayer(wafer,indx.first); +} + +int HGCalDDDConstants::wafers() const { + + int wafer(0); + for (unsigned int i = 0; idepth_[i]; + wafer += modules(lay, true); } - for (unsigned int i=0; i& pos) const { + int num(0); + const double tol(0.00001); + double cellY = 2.0*cellR*tan30deg_; + for (unsigned int k=0; k 0) { - trform_.back().h3v /= nz; - } - } - } - } -#ifdef DebugLog - std::cout << "Obtained " << trform_.size() << " transformation matrices" - << std::endl; - for (unsigned int k=0; k illegal"; - throw cms::Exception("DDException") << "HGCalDDDConstants: cannot get array " << str; - } - } else { - if (nval < 1 && nmin == 0) { - edm::LogError("HGCalGeom") << "HGCalDDDConstants : # of " << str - << " bins " << nval << " < 2 ==> illegal" - << " (nmin=" << nmin << ")"; - throw cms::Exception("DDException") << "HGCalDDDConstants: cannot get array " << str; - } - } - nmin = nval; - return fvec; +void HGCalDDDConstants::getParameterSquare(int lay, int subSec, bool reco, + float& h, float& bl, float& tl, + float& alpha) const { + if (reco) { + h = hgpar_->moduler_[lay].h; + bl = hgpar_->moduler_[lay].bl; + tl = hgpar_->moduler_[lay].tl; + alpha= hgpar_->moduler_[lay].alpha; + if ((subSec>0 && alpha<0) || (subSec<=0 && alpha>0)) alpha = -alpha; } else { - if (nmin >= 0) { - edm::LogError("HGCalGeom") << "HGCalDDDConstants: cannot get array " - << str; - throw cms::Exception("DDException") << "HGCalDDDConstants: cannot get array " << str; - } - std::vector fvec; - nmin = 0; - return fvec; + h = hgpar_->modules_[lay].h; + bl = hgpar_->modules_[lay].bl; + tl = hgpar_->modules_[lay].tl; + alpha= hgpar_->modules_[lay].alpha; } } -std::pair HGCalDDDConstants::getIndex(int lay, bool reco) const { +bool HGCalDDDConstants::waferInLayer(int wafer, int lay) const { - if (lay<1 || lay>(int)(layerIndex.size())) return std::pair(-1,0); - if (reco && lay>(int)(depthIndex.size())) return std::pair(-1,0); - int i = (reco ? depthIndex[lay-1] : layerIndex[lay-1]); - float cell = (reco ? moduler_[i].cellSize : modules_[i].cellSize); - return std::pair(i,cell); + double rr = 2*rmax_*tan30deg_; + double rpos = hgpar_->waferPos_[wafer].perp(); + bool in = (rpos-rr >= hgpar_->rMinLayHex_[lay] && + rpos+rr <= hgpar_->rMaxLayHex_[lay]); + return in; } #include "FWCore/Utilities/interface/typelookup.h" diff --git a/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc b/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc new file mode 100644 index 0000000000000..56e8c51cccea3 --- /dev/null +++ b/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc @@ -0,0 +1,752 @@ +#include "Geometry/HGCalCommonData/interface/HGCalGeomParameters.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "FWCore/Utilities/interface/Exception.h" + +#include "DetectorDescription/Base/interface/DDutils.h" +#include "DetectorDescription/Core/interface/DDValue.h" +#include "DetectorDescription/Core/interface/DDFilter.h" +#include "DetectorDescription/Core/interface/DDSolid.h" +#include "DetectorDescription/Core/interface/DDConstant.h" +#include "DetectorDescription/Core/interface/DDVectorGetter.h" +#include "DetectorDescription/Core/interface/DDFilteredView.h" +#include "DetectorDescription/RegressionTest/interface/DDErrorDetection.h" +#include "CondFormats/GeometryObjects/interface/HGCalParameters.h" +#include "CLHEP/Units/GlobalPhysicalConstants.h" +#include "CLHEP/Units/GlobalSystemOfUnits.h" + +//#define DebugLog + +const double k_ScaleFromDDD = 0.1; + +HGCalGeomParameters::HGCalGeomParameters() { +#ifdef DebugLog + std::cout << "HGCalGeomParameters::HGCalGeomParameters() constructor" << std::endl; +#endif +} + +HGCalGeomParameters::~HGCalGeomParameters() { +#ifdef DebugLog + std::cout << "HGCalGeomParameters::destructed!!!" << std::endl; +#endif +} + +void HGCalGeomParameters::loadGeometrySquare(const DDFilteredView& _fv, + HGCalParameters& php, + const std::string & sdTag) { + + DDFilteredView fv = _fv; + bool dodet(true), first(true); + int zpFirst(0); + std::vector trforms; + + while (dodet) { + const DDSolid & sol = fv.logicalPart().solid(); + std::string name = sol.name(); + int isd = (name.find(sdTag) == std::string::npos) ? -1 : 1; + if (isd > 0) { + std::vector copy = fv.copyNumbers(); + int nsiz = (int)(copy.size()); + int lay = (nsiz > 0) ? copy[nsiz-1] : -1; + int sec = (nsiz > 1) ? copy[nsiz-2] : -1; + int zp = (nsiz > 3) ? copy[nsiz-4] : -1; + if (zp !=1 ) zp = -1; + if (first) {first = false; zpFirst = zp;} + const DDTrap & trp = static_cast(sol); + HGCalParameters::hgtrap mytr; + mytr.lay = lay; mytr.bl = trp.x1(); mytr.tl = trp.x2(); + mytr.h = 0.5*(trp.y1()+trp.y2()); mytr.dz = trp.halfZ(); + mytr.alpha = trp.alpha1(); mytr.cellSize = 0; + int subs = (trp.alpha1()>0 ? 1 : 0); + if (std::find(php.layer_.begin(),php.layer_.end(),lay) == + php.layer_.end()) { + for (unsigned int k=0; k 0) { + php.trform_.back().h3v /= nz; + } + } + } + } +#ifdef DebugLog + std::cout << "Obtained " << php.trform_.size() << " transformation matrices" + << std::endl; + for (unsigned int k=0; k layers; + std::vector trforms; + + while (dodet) { + const DDSolid & sol = fv.logicalPart().solid(); + std::string name = sol.name(); + bool isd = (name.find(sdTag1) != std::string::npos); + // Layers first + if (isd) { + std::vector copy = fv.copyNumbers(); + int nsiz = (int)(copy.size()); + int lay = (nsiz > 0) ? copy[nsiz-1] : 0; + int zp = (nsiz > 2) ? copy[nsiz-3] : -1; + if (zp !=1 ) zp = -1; + if (lay == 0) { + edm::LogError("HGCalGeom") << "Funny layer # " << lay << " zp " + << zp << " in " << nsiz << " components"; + throw cms::Exception("DDException") << "Funny layer # " << lay; + } else { + if (std::find(php.layer_.begin(),php.layer_.end(),lay) == + php.layer_.end()) php.layer_.push_back(lay); + std::map::iterator itr = layers.find(lay); + if (itr == layers.end()) { + const DDTubs & tube = static_cast(sol); + double rin = k_ScaleFromDDD*tube.rIn(); + double rout= k_ScaleFromDDD*tube.rOut(); + double zp = k_ScaleFromDDD*fv.translation().Z(); + HGCalGeomParameters::layerParameters laypar(rin,rout,zp); + layers[lay] = laypar; + } + DD3Vector x, y, z; + fv.rotation().GetComponents( x, y, z ) ; + const CLHEP::HepRep3x3 rotation ( x.X(), y.X(), z.X(), + x.Y(), y.Y(), z.Y(), + x.Z(), y.Z(), z.Z() ); + const CLHEP::HepRotation hr ( rotation ); + double xx = k_ScaleFromDDD*fv.translation().X(); + if (std::abs(xx) < 0.001) xx = 0; + double yy = k_ScaleFromDDD*fv.translation().Y(); + if (std::abs(yy) < 0.001) yy = 0; + const CLHEP::Hep3Vector h3v ( xx, yy, fv.translation().Z() ); + HGCalParameters::hgtrform mytrf; + mytrf.zp = zp; + mytrf.lay = lay; + mytrf.sec = 0; + mytrf.subsec= 0; + mytrf.h3v = h3v; + mytrf.hr = hr; + mytrf.used = false; + trforms.push_back(mytrf); + } + } + dodet = fv.next(); + } + + // Then wafers + std::map wafers; + std::string attribute = "Volume"; + DDValue val1(attribute, sdTag2, 0.0); + DDSpecificsFilter filter1; + filter1.setCriteria(val1, DDCompOp::equals); + DDFilteredView fv1(*cpv); + fv1.addFilter(filter1); + bool ok = fv1.firstChild(); + double rmax(0); + if (!ok) { + edm::LogError("HGCalGeom") << " Attribute " << val1 + << " not found but needed."; + throw cms::Exception("DDException") << "Attribute " << val1 + << " not found but needed."; + } else { + dodet = true; + std::vector names; + while (dodet) { + const DDSolid & sol = fv1.logicalPart().solid(); + std::string name = fv1.logicalPart().name(); + bool isd = (name.find(sdTag2) != std::string::npos); + if (isd) { + std::vector copy = fv1.copyNumbers(); + int nsiz = (int)(copy.size()); + int wafer= (nsiz > 0) ? copy[nsiz-1] : -1; + if (wafer == 0) { + edm::LogError("HGCalGeom") << "Funny wafer # " << wafer << " in " + << nsiz << " components"; + throw cms::Exception("DDException") << "Funny wafer # " << wafer; + } else { + std::map::iterator itr = wafers.find(wafer); + if (itr == wafers.end()) { + double xx = k_ScaleFromDDD*fv1.translation().X(); + if (std::abs(xx) < 0.001) xx = 0; + double yy = k_ScaleFromDDD*fv1.translation().Y(); + if (std::abs(yy) < 0.001) yy = 0; + wafers[wafer] = GlobalPoint(xx,yy,k_ScaleFromDDD*fv1.translation().Z()); + std::vector::iterator its = std::find(names.begin(),names.end(),name); + if (its == names.end()) { + const DDPolyhedra & polyhedra = static_cast(sol); + std::vector zv = polyhedra.zVec(); + std::vector rv = polyhedra.rMaxVec(); + php.waferR_ = rv[0]/std::cos(30.0*CLHEP::deg); + rmax = k_ScaleFromDDD*rv[0]; + double dz = 0.5*(zv[1]-zv[0]); + HGCalParameters::hgtrap mytr; + mytr.lay = 1; mytr.bl = php.waferR_; + mytr.tl = php.waferR_; mytr.h = php.waferR_; + mytr.dz = dz; mytr.alpha = 0.0; + mytr.cellSize = waferSize_; + php.modules_.push_back(mytr); + names.push_back(name); + } + } + } + } + dodet = fv1.next(); + } + } + + // Finally the cells + std::map wafertype; + std::map cellsf, cellsc; + double ymax = 2.0*rmax*tan(30.0*CLHEP::deg); + DDValue val2(attribute, sdTag3, 0.0); + DDSpecificsFilter filter2; + filter2.setCriteria(val2, DDCompOp::equals); + DDFilteredView fv2(*cpv); + fv2.addFilter(filter2); + ok = fv2.firstChild(); + if (!ok) { + edm::LogError("HGCalGeom") << " Attribute " << val2 + << " not found but needed."; + throw cms::Exception("DDException") << "Attribute " << val2 + << " not found but needed."; + } else { + dodet = true; + while (dodet) { + const DDSolid & sol = fv2.logicalPart().solid(); + std::string name = sol.name(); + bool isd = (name.find(sdTag3) != std::string::npos); + if (isd) { + std::vector copy = fv2.copyNumbers(); + int nsiz = (int)(copy.size()); + int cellx= (nsiz > 0) ? copy[nsiz-1] : 0; + int wafer= (nsiz > 1) ? copy[nsiz-2] : 0; + int cell = cellx%1000; + int type = cellx/1000; + if (type != 1 && type != 2) { + edm::LogError("HGCalGeom") << "Funny cell # " << cell << " type " + << type << " in " << nsiz << " components"; + throw cms::Exception("DDException") << "Funny cell # " << cell; + } else { + std::map::iterator ktr = wafertype.find(wafer); + if (ktr == wafertype.end()) wafertype[wafer] = type; + bool newc(false); + std::map::iterator itr; + if (type == 1) { + itr = cellsf.find(cell); + newc= (itr == cellsf.end()); + } else { + itr = cellsc.find(cell); + newc= (itr == cellsc.end()); + } + if (newc) { + bool half = (name.find("Half") != std::string::npos); + double xx = k_ScaleFromDDD*fv2.translation().X(); + double yy = k_ScaleFromDDD*fv2.translation().Y(); + HGCalGeomParameters::cellParameters cp(half,GlobalPoint(xx,yy,0)); + if (type == 1) { + cellsf[cell] = cp; + } else { + cellsc[cell] = cp; + } + } + } + } + dodet = fv2.next(); + } + } + + if ((cellsf.size()==0) || (cellsc.size()==0) || (wafers.size()==0) || + (layers.size()==0)) { + edm::LogError("HGCalGeom") << "HGCalGeomParameters : number of cells " + << cellsf.size() << ":" << cellsc.size() + << " wafers " << wafers.size() << " layers " + << layers.size() << " illegal"; + throw cms::Exception("DDException") + << "HGCalGeomParameters: mismatch between geometry and specpar: cells " + << cellsf.size() << ":" << cellsc.size() << " wafers " << wafers.size() + << " layers " << layers.size(); + } + + for (unsigned int i=0; i::iterator itr = layers.begin(); + itr != layers.end(); ++itr) { + if (itr->first == (int)(i+1)) { + php.layerIndex_.push_back(itr->first); + php.rMinLayHex_.push_back(itr->second.rmin); + php.rMaxLayHex_.push_back(itr->second.rmax); + php.zLayerHex_.push_back(itr->second.zpos); + break; + } + } + } + for (unsigned int i=0; i 0) { + php.trform_.back().h3v /= nz; + } + } + } + } + + double rmin = k_ScaleFromDDD*php.waferR_; + for (std::map::iterator itr = wafers.begin(); + itr != wafers.end(); ++itr) { + php.waferCopy_.push_back(itr->first); + php.waferPos_.push_back(itr->second); + std::map::iterator ktr = wafertype.find(itr->first); + int typet = (ktr == wafertype.end()) ? 0 : (ktr->second); + php.waferTypeT_.push_back(typet); + double r = (itr->second).perp(); + int type(3); + for (int k=1; k<4; ++k) { + if ((r+rmin)<=php.boundR_[k]) { + type = k; break; + } + } + php.waferTypeL_.push_back(type); + } + php.nSectors_ = (int)(php.waferCopy_.size()); + + std::map::const_iterator itrf = wafers.end(); + for (unsigned int i=0; i::iterator itr = cellsf.find(i); + if (itr == cellsf.end()) { + edm::LogError("HGCalGeom") << "HGCalGeomParameters: missing info for" + << " fine cell number " << i; + throw cms::Exception("DDException") + << "HGCalGeomParameters: missing info for fine cell number " << i; + } else { + double xx = (itr->second).xyz.x(); + double yy = (itr->second).xyz.y(); + std::pair xy = cellPosition(wafers,itrf,i,rmax,ymax,xx,yy,cellsf.size()); + if ((itr->second).half) { + if (xy.first > 0) xy.first -= 0.001; + else xy.first += 0.001; + } + php.cellFine_.push_back(GlobalPoint(xy.first,xy.second,0)); + } + } + itrf = wafers.end(); + for (unsigned int i=0; i::iterator itr = cellsc.find(i); + if (itr == cellsc.end()) { + edm::LogError("HGCalGeom") << "HGCalGeomParameters: missing info for" + << " coarse cell number " << i; + throw cms::Exception("DDException") + << "HGCalGeomParameters: missing info for coarse cell number " << i; + } else { + double xx = (itr->second).xyz.x(); + double yy = (itr->second).xyz.y(); + std::pair xy = cellPosition(wafers,itrf,i,rmax,ymax,xx,yy,cellsc.size()); + if ((itr->second).half) { + if (xy.first > 0) xy.first -= 0.001; + else xy.first += 0.001; + } + php.cellCoarse_.push_back(GlobalPoint(xy.first,xy.second,0)); + } + } + int depth(0); + for (unsigned int i=0; i dummy = getDDDArray("WaferSize",sv,nmin); + waferSize_ = dummy[0]; + } +#ifdef DebugLog + std::cout << "HGCalGeomParameters: Wafer Size: " << waferSize_ << std::endl; +#endif + + //Cell size + DDValue val2(attribute, sdTag2, 0.0); + DDSpecificsFilter filter2; + filter2.setCriteria(val2, DDCompOp::equals); + DDFilteredView fv2(*cpv); + fv2.addFilter(filter2); + if (fv2.firstChild()) { + DDsvalues_type sv(fv2.mergedSpecifics()); + int nmin(0); + php.cellSize_ = getDDDArray("CellSize",sv,nmin); + } +#ifdef DebugLog + std::cout << "HGCalGeomParameters: " << php.cellSize_.size() + << " cells of sizes:"; + for (unsigned int k=0; k HGCalGeomParameters::getDDDArray(const std::string & str, + const DDsvalues_type & sv, + int & nmin) { + DDValue value(str); + if (DDfetch(&sv,value)) { + const std::vector & fvec = value.doubles(); + int nval = fvec.size(); + if (nmin > 0) { + if (nval < nmin) { + edm::LogError("HGCalGeom") << "HGCalGeomParameters : # of " << str + << " bins " << nval << " < " << nmin + << " ==> illegal"; + throw cms::Exception("DDException") << "HGCalGeomParameters: cannot get array " << str; + } + } else { + if (nval < 1 && nmin == 0) { + edm::LogError("HGCalGeom") << "HGCalGeomParameters : # of " << str + << " bins " << nval << " < 2 ==> illegal" + << " (nmin=" << nmin << ")"; + throw cms::Exception("DDException") << "HGCalGeomParameters: cannot get array " << str; + } + } + nmin = nval; + return fvec; + } else { + if (nmin >= 0) { + edm::LogError("HGCalGeom") << "HGCalGeomParameters: cannot get array " + << str; + throw cms::Exception("DDException") << "HGCalGeomParameters: cannot get array " << str; + } + std::vector fvec; + nmin = 0; + return fvec; + } +} + +std::pair +HGCalGeomParameters::cellPosition(const std::map& wafers, + std::map::const_iterator& itrf, + unsigned int num, double rmax, double ymax, + double xx, double yy, unsigned int ncells) { + + if (itrf == wafers.end()) { + for (std::map::const_iterator itr = wafers.begin(); + itr != wafers.end(); ++itr) { + double dx = std::abs(xx - itr->second.x()); + double dy = std::abs(yy - itr->second.y()); + if (dx <= (rmax+0.0001) && dy <= ymax) { + double xmax = (dy <= 0.5*ymax) ? rmax : (rmax-(dy-0.5*ymax)/tan(30.0*CLHEP::deg)); + if ((dx <= (xmax+0.0001)) && ((yy-itr->second.y()>0) || (num<=ncells/2+8))) { + itrf = itr; + break; + } + } + } + } + double dx(0), dy(0); + if (itrf != wafers.end()) { + dx = (xx - itrf->second.x()); + if (std::abs(dx) < 0.001) dx = 0; + dy = (yy - itrf->second.y()); + if (std::abs(dy) < 0.001) dy = 0; + } + return std::pair(dx,dy); +} diff --git a/Geometry/HGCalCommonData/src/HGCalGeometryMode.cc b/Geometry/HGCalCommonData/src/HGCalGeometryMode.cc new file mode 100644 index 0000000000000..df7e44566662d --- /dev/null +++ b/Geometry/HGCalCommonData/src/HGCalGeometryMode.cc @@ -0,0 +1,8 @@ +#include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h" + +template<> +HGCalStringToEnumParser::HGCalStringToEnumParser() { + enumMap["HGCalGeometryMode::Square"] = HGCalGeometryMode::Square; + enumMap["HGCalGeometryMode::Hexagon"] = HGCalGeometryMode::Hexagon; +} + diff --git a/Geometry/HGCalCommonData/src/HGCalParametersFromDD.cc b/Geometry/HGCalCommonData/src/HGCalParametersFromDD.cc new file mode 100644 index 0000000000000..b860be7ba725c --- /dev/null +++ b/Geometry/HGCalCommonData/src/HGCalParametersFromDD.cc @@ -0,0 +1,90 @@ +#include "Geometry/HGCalCommonData/interface/HGCalParametersFromDD.h" +#include "Geometry/HGCalCommonData/interface/HGCalGeomParameters.h" +#include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h" +#include "CondFormats/GeometryObjects/interface/HGCalParameters.h" +#include "DetectorDescription/Core/interface/DDCompactView.h" +#include "DetectorDescription/Core/interface/DDFilteredView.h" +#include "DetectorDescription/Core/interface/DDVectorGetter.h" +#include "DetectorDescription/Base/interface/DDutils.h" + +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "CLHEP/Units/GlobalSystemOfUnits.h" +#include +#include + +//#define DebugLog + +namespace { + int getGeometryMode(const char* s, const DDsvalues_type & sv) { + DDValue val(s); + if (DDfetch(&sv, val)) { + const std::vector & fvec = val.strings(); + if (fvec.size() == 0) { + throw cms::Exception("HGCalGeom") << "Failed to get " << s << " tag."; + } + + int result(-1); + HGCalStringToEnumParser eparser; + HGCalGeometryMode::GeometryMode mode = (HGCalGeometryMode::GeometryMode) eparser.parseString(fvec[0]); + result = (int)(mode); + return result; + } else { + throw cms::Exception("HGCalGeom") << "Failed to get "<< s << " tag."; + } + } +} + +bool HGCalParametersFromDD::build(const DDCompactView* cpv, + HGCalParameters& php, + const std::string& name, + const std::string& namew, + const std::string& namec) { + +#ifdef DebugLog + std::cout << "HGCalParametersFromDD::build called with names " << name << ":" + << namew << ":" << namec << std::endl; +#endif + + //Special parameters at simulation level + std::string attribute = "Volume"; + std::string value = name; + DDValue val(attribute, value, 0.0); + DDSpecificsFilter filter; + filter.setCriteria(val, DDCompOp::equals); + DDFilteredView fv(*cpv); + fv.addFilter(filter); + bool ok = fv.firstChild(); + + if (ok) { + DDsvalues_type sv(fv.mergedSpecifics()); + php.mode_ = getGeometryMode("GeometryMode", sv); + HGCalGeomParameters *geom = new HGCalGeomParameters(); + if (php.mode_ == HGCalGeometryMode::Square) { + //Load the SpecPars + geom->loadSpecParsSquare(fv, php); + //Load the Geometry parameters + geom->loadGeometrySquare(fv, php, name); + } else if (php.mode_ == HGCalGeometryMode::Hexagon) { + //Load the SpecPars + geom->loadSpecParsHexagon(fv, php, cpv, namew, namec); + //Load the Geometry parameters + geom->loadGeometryHexagon(fv, php, name, cpv, namew, namec); + } else { + edm::LogError("HGCalGeom") << "Unknown Geometry type " << php.mode_ + << " for HGCal " << name << ":" << namew + << ":" << namec; + throw cms::Exception("DDException") + << "Unknown Geometry type " << php.mode_ << " for HGCal " << name + << ":" << namew << ":" << namec; + } + } else { + edm::LogError("HGCalGeom") << " Attribute " << val + << " not found but needed."; + throw cms::Exception("DDException") << "Attribute " << val + << " not found but needed."; + } + + edm::LogInfo("HGCalGeom") << "Return from HGCalParametersFromDD::build with " + << ok; + return ok; +} diff --git a/Geometry/HGCalCommonData/test/HGCGeometryTester.cc b/Geometry/HGCalCommonData/test/HGCGeometryTester.cc index 9de91e01fec90..9b1922ce02565 100644 --- a/Geometry/HGCalCommonData/test/HGCGeometryTester.cc +++ b/Geometry/HGCalCommonData/test/HGCGeometryTester.cc @@ -50,9 +50,13 @@ class HGCGeometryTester : public edm::one::EDAnalyzer<> { void beginJob() override {} void analyze(edm::Event const& iEvent, edm::EventSetup const&) override; void endJob() override {} +private: + bool square; }; -HGCGeometryTester::HGCGeometryTester(const edm::ParameterSet& ) {} +HGCGeometryTester::HGCGeometryTester(const edm::ParameterSet& iC) { + square = iC.getUntrackedParameter("SquareType",true); +} HGCGeometryTester::~HGCGeometryTester() {} @@ -79,10 +83,17 @@ void HGCGeometryTester::analyze( const edm::Event& iEvent, if (svPars.find(name) == svPars.end()) { //print half height and widths for the trapezoid std::vector solidPar=eview.logicalPart().solid().parameters(); - svPars[name] = std::pair(solidPar[3], - 0.5*(solidPar[4]+solidPar[5])); - std::cout << name << " Layer " << layer << " " << solidPar[3] - << " " << solidPar[4] << " " << solidPar[5] << std::endl; + if (square) { + svPars[name] = std::pair(solidPar[3], + 0.5*(solidPar[4]+solidPar[5])); + std::cout << name << " Layer " << layer << " " << solidPar[3] + << " " << solidPar[4] << " " << solidPar[5] << std::endl; + } else { + svPars[name] = std::pair(solidPar[0], + 0.5*(solidPar[2]-solidPar[1])); + std::cout << name << " Layer " << layer << " " << solidPar[0] + << " " << solidPar[1] << " " << solidPar[2] << std::endl; + } } } }while(eview.next() ); diff --git a/Geometry/HGCalCommonData/test/HGCalNumberingTester.cc b/Geometry/HGCalCommonData/test/HGCalNumberingTester.cc index 9a79b19dfa255..e2d0513e07363 100644 --- a/Geometry/HGCalCommonData/test/HGCalNumberingTester.cc +++ b/Geometry/HGCalCommonData/test/HGCalNumberingTester.cc @@ -13,7 +13,6 @@ // // Original Author: Sunanda Banerjee // Created: Mon 2014/03/21 -// $Id: HGCalNumberingTester.cc,v 1.0 2014/032/21 14:06:07 sunanda Exp $ // // @@ -22,6 +21,7 @@ #include #include #include +#include // user include files #include "FWCore/Framework/interface/Frameworkfwd.h" @@ -40,10 +40,10 @@ #include "Geometry/Records/interface/IdealGeometryRecord.h" #include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h" +#include "CLHEP/Units/GlobalSystemOfUnits.h" #include "CoralBase/Exception.h" -class HGCalNumberingTester : public edm::one::EDAnalyzer<> -{ +class HGCalNumberingTester : public edm::one::EDAnalyzer<> { public: explicit HGCalNumberingTester( const edm::ParameterSet& ); ~HGCalNumberingTester(); @@ -51,9 +51,21 @@ class HGCalNumberingTester : public edm::one::EDAnalyzer<> void beginJob() override {} void analyze(edm::Event const& iEvent, edm::EventSetup const&) override; void endJob() override {} +private: + std::string nameSense_, nameDetector_; + double position_; + int increment_; }; -HGCalNumberingTester::HGCalNumberingTester(const edm::ParameterSet& ) {} +HGCalNumberingTester::HGCalNumberingTester(const edm::ParameterSet& iC) { + nameSense_ = iC.getParameter("NameSense"); + nameDetector_ = iC.getParameter("NameDevice"); + position_ = iC.getParameter("LocalPosition")*CLHEP::mm; + increment_ = iC.getParameter("Increment"); + std::cout << "Test numbering for " << nameDetector_ <<" using constants of " + << nameSense_ << " at local position of " << position_/CLHEP::mm + << " mm for every " << increment_ << " layers" << std::endl; +} HGCalNumberingTester::~HGCalNumberingTester() {} @@ -62,109 +74,44 @@ void HGCalNumberingTester::analyze( const edm::Event& iEvent, const edm::EventSe edm::ESHandle pHGNDC; - iSetup.get().get("HGCalEESensitive",pHGNDC); - const HGCalDDDConstants hgeedc(*pHGNDC); - std::cout << "EE Layers = " << hgeedc.layers(false) << " Sectors = " - << hgeedc.sectors() << std::endl; + iSetup.get().get(nameSense_,pHGNDC); + const HGCalDDDConstants hgdc(*pHGNDC); + std::cout << nameDetector_ << " Layers = " << hgdc.layers(false) + << " Sectors = " << hgdc.sectors() << std::endl; std::pair kxy, lxy; std::pair xy; - float localx(5.0), localy(5.0); - for (unsigned int i=0; i ncells = hgeedc.numberCells(i+1,false); - std::cout << "Layer " << i+1 << " with " << ncells.size() << " rows\n"; - int ntot(0); - for (unsigned int k=0; k().get("HGCalHESiliconSensitive",pHGNDC); - const HGCalDDDConstants hghesidc(*pHGNDC); - std::cout << "HE Silicon Layers = " << hghesidc.layers(false) - << " Sectors = " << hghesidc.sectors() << std::endl; - for (unsigned int i=0; i ncells = hghesidc.numberCells(i+1,false); - std::cout << "Layer " << i+1 << " with " << ncells.size() << " rows\n"; - int ntot(0); - for (unsigned int k=0; k().get("HGCalHEScintillatorSensitive",pHGNDC); - const HGCalDDDConstants hghescdc(*pHGNDC); - std::cout << "HE Scintillator Layers = " << hghescdc.layers(false) - << " Sectors = " << hghescdc.sectors() << std::endl; - std::vector::const_iterator itr = hghescdc.getFirstModule(false); - int subsec = ((itr->alpha) > 0) ? 1 : 0; - for (unsigned int i=0; i::const_iterator itr = hgdc.getFirstModule(false); + bool halfCell = ((itr->alpha) > 0); + int subsec = (halfCell) ? 1 : 0; + float localx(position_), localy(position_); + for (unsigned int i=0; i ncells = hghescdc.numberCells(i+1,false); + << ", " << subsec << "), assignCell o/p (" << kxy.first << ", " + << kxy.second << ") locateCell o/p (" << xy.first << ", " + << xy.second << ")," << " final (" << lxy.first << ", " + << lxy.second << ")" << std::endl; + std::vector ncells = hgdc.numberCells(i+1,false); std::cout << "Layer " << i+1 << " with " << ncells.size() << " rows\n"; int ntot(0); for (unsigned int k=0; k ParmVec; + HGCalGeometryLoader (); ~HGCalGeometryLoader (); HGCalGeometry* build(const HGCalTopology& ); + +private: + void buildGeom(const ParmVec&, const HepGeom::Transform3D&, const DetId&, + HGCalGeometry*); + }; #endif diff --git a/Geometry/HGCalGeometry/plugins/HGCalGeometryESProducer.cc b/Geometry/HGCalGeometry/plugins/HGCalGeometryESProducer.cc index b81ed3224278b..d3ed22c55d9f5 100644 --- a/Geometry/HGCalGeometry/plugins/HGCalGeometryESProducer.cc +++ b/Geometry/HGCalGeometry/plugins/HGCalGeometryESProducer.cc @@ -30,7 +30,7 @@ #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h" #include "Geometry/HGCalGeometry/interface/HGCalGeometryLoader.h" -//#define DebugLog +#define DebugLog // // class decleration diff --git a/Geometry/HGCalGeometry/src/HGCalGeometry.cc b/Geometry/HGCalGeometry/src/HGCalGeometry.cc index 89b3825de32ba..6f5e5d2f9b4e1 100644 --- a/Geometry/HGCalGeometry/src/HGCalGeometry.cc +++ b/Geometry/HGCalGeometry/src/HGCalGeometry.cc @@ -26,8 +26,8 @@ HGCalGeometry::HGCalGeometry( const HGCalTopology& topology_ ) m_cellVec( topology_.totalGeomModules()), m_validGeomIds( topology_.totalGeomModules()), m_halfType( topology_.detectorType()), - m_subdet( topology_.subDetector()) -{ + m_subdet( topology_.subDetector()) { + m_validIds.reserve( topology().totalModules()); #ifdef DebugLog std::cout << "Expected total # of Geometry Modules " @@ -54,19 +54,32 @@ void HGCalGeometry::newCell( const GlobalPoint& f1 , const GlobalPoint& f3 , const CCGFloat* parm , const DetId& detId ) { + + DetId geomId; + int cells; + HGCalTopology::DecodedDetId id = topology().decode(detId); + if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) { + geomId = (detId.subdetId() == HGCEE ? + (DetId)(HGCEEDetId(detId).geometryCell()) : + (DetId)(HGCHEDetId(detId).geometryCell())); + cells = topology().dddConstants().maxCells(id.iLay,true); + } else { + geomId = (DetId)(HGCalDetId(detId).geometryCell()); + cells = topology().dddConstants().numberCellsHexagon(id.iSec); + } const uint32_t cellIndex (topology().detId2denseGeomId(detId)); - DetId geomId = (detId.subdetId() == HGCEE ? - (DetId)(HGCEEDetId(detId).geometryCell()) : - (DetId)(HGCHEDetId(detId).geometryCell())); + m_cellVec.at( cellIndex ) = FlatTrd( cornersMgr(), f1, f2, f3, parm ) ; m_validGeomIds.at( cellIndex ) = geomId ; - HGCalTopology::DecodedDetId id = topology().decode(detId); - int cells = topology().dddConstants().maxCells(id.iLay,true); +#ifdef DebugLog + unsigned int nOld = m_validIds.size(); +#endif for (int cell = 0; cell < cells; ++cell) { id.iCell = cell; m_validIds.push_back( topology().encode(id)); - if (!m_halfType) { + if ((topology().dddConstants().geomMode() == HGCalGeometryMode::Square) && + (!m_halfType)) { id.iSubSec = -id.iSubSec; m_validIds.push_back( topology().encode(id)); id.iSubSec = -id.iSubSec; @@ -79,8 +92,14 @@ void HGCalGeometry::newCell( const GlobalPoint& f1 , << " back:" << f2.x() << '/' << f2.y() << '/' << f2.z() << " eta|phi " << m_cellVec[cellIndex].etaPos() << ":" << m_cellVec[cellIndex].phiPos() << " id:"; - if (m_subdet == HGCEE) std::cout << HGCEEDetId(detId); - else std::cout << HGCHEDetId(detId); + if (topology().dddConstants().geomMode() != HGCalGeometryMode::Square) { + std::cout << HGCalDetId(detId); + } else if (m_subdet == HGCEE) { + std::cout << HGCEEDetId(detId); + } else { + std::cout << HGCHEDetId(detId); + } + unsigned int nNew = m_validIds.size(); std::cout << " with valid DetId from " << nOld << " to " << nNew << std::endl; std::cout << "Cell[" << cellIndex << "] " << std::hex << geomId.rawId() @@ -91,10 +110,16 @@ void HGCalGeometry::newCell( const GlobalPoint& f1 , const CaloCellGeometry* HGCalGeometry::getGeometry(const DetId& id) const { if (id == DetId()) return 0; // nothing to get - DetId geoId = (id.subdetId() == HGCEE ? - (DetId)(HGCEEDetId(id).geometryCell()) : - (DetId)(HGCHEDetId(id).geometryCell())); + DetId geoId; + if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) { + geoId = (id.subdetId() == HGCEE ? + (DetId)(HGCEEDetId(id).geometryCell()) : + (DetId)(HGCHEDetId(id).geometryCell())); + } else { + geoId = (DetId)(HGCalDetId(id).geometryCell()); + } const uint32_t cellIndex (topology().detId2denseGeomId(geoId)); + std::cout << HGCalDetId(geoId) << " cellIndex " << cellIndex << std::endl; /* if (cellIndex < m_cellVec.size()) { HGCalTopology::DecodedDetId id_ = topology().decode(id); @@ -108,12 +133,17 @@ const CaloCellGeometry* HGCalGeometry::getGeometry(const DetId& id) const { } -GlobalPoint HGCalGeometry::getPosition( const DetId& id ) const { +GlobalPoint HGCalGeometry::getPosition(const DetId& id) const { unsigned int cellIndex = indexFor(id); if (cellIndex < m_cellVec.size()) { HGCalTopology::DecodedDetId id_ = topology().decode(id); - std::pair xy = topology().dddConstants().locateCell(id_.iCell,id_.iLay,id_.iSubSec,true); + std::pair xy; + if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) { + xy = topology().dddConstants().locateCell(id_.iCell,id_.iLay,id_.iSubSec,true); + } else { + xy = topology().dddConstants().locateCellHex(id_.iCell,id_.iSec,true); + } const HepGeom::Point3D lcoord(xy.first,xy.second,0); #ifdef DebugLog std::cout << "getPosition:: index " << cellIndex << " Local " << xy.first @@ -126,13 +156,18 @@ GlobalPoint HGCalGeometry::getPosition( const DetId& id ) const { return GlobalPoint(); } -HGCalGeometry::CornersVec HGCalGeometry::getCorners( const DetId& id ) const { +HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& id) const { HGCalGeometry::CornersVec co (8, GlobalPoint(0,0,0)); unsigned int cellIndex = indexFor(id); if (cellIndex < m_cellVec.size()) { HGCalTopology::DecodedDetId id_ = topology().decode(id); - std::pair xy = topology().dddConstants().locateCell(id_.iCell,id_.iLay,id_.iSubSec,true); + std::pair xy; + if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) { + xy = topology().dddConstants().locateCell(id_.iCell,id_.iLay,id_.iSubSec,true); + } else { + xy = topology().dddConstants().locateCellHex(id_.iCell,id_.iSec,true); + } float dz = m_cellVec[cellIndex].param()[0]; float dx = 0.5*m_cellVec[cellIndex].param()[11]; static const int signx[] = {-1,-1,1,1,-1,-1,1,1}; @@ -146,16 +181,28 @@ HGCalGeometry::CornersVec HGCalGeometry::getCorners( const DetId& id ) const { return co; } -DetId HGCalGeometry::getClosestCell( const GlobalPoint& r ) const { +DetId HGCalGeometry::getClosestCell(const GlobalPoint& r) const { unsigned int cellIndex = getClosestCellIndex(r); if (cellIndex < m_cellVec.size()) { - const HepGeom::Point3D local = m_cellVec[cellIndex].getLocal(r); HGCalTopology::DecodedDetId id_ = topology().decode(m_validGeomIds[cellIndex]); + HepGeom::Point3D local; + if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) { + local = m_cellVec[cellIndex].getLocal(r); + } else { + if (r.z() > 0) local = HepGeom::Point3D(r.x(),r.y(),0); + else local = HepGeom::Point3D(-r.x(),r.y(),0); + } std::pair kxy = topology().dddConstants().assignCell(local.x(),local.y(),id_.iLay, id_.iSubSec,true); id_.iCell = kxy.second; - id_.iSubSec = kxy.first; + if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) { + id_.iSubSec = kxy.first; + } else { + id_.iSec = kxy.first; + id_.iSubSec = topology().dddConstants().waferTypeT(kxy.first); + if (id_.iSubSec != 1) id_.iSubSec = -1; + } #ifdef DebugLog std::cout << "getClosestCell: local " << local << " Id " << id_.zside << ":" << id_.iLay << ":" << id_.iSec << ":" << id_.iSubSec @@ -163,7 +210,7 @@ DetId HGCalGeometry::getClosestCell( const GlobalPoint& r ) const { #endif //check if returned cell is valid - if(id_.iCell>=0) return topology().encode(id_); + if (id_.iCell>=0) return topology().encode(id_); } //if not valid or out of bounds return a null DetId @@ -184,9 +231,14 @@ std::string HGCalGeometry::cellElement() const { unsigned int HGCalGeometry::indexFor(const DetId& id) const { unsigned int cellIndex = m_cellVec.size(); if (id != DetId()) { - DetId geoId = (id.subdetId() == HGCEE ? - (DetId)(HGCEEDetId(id).geometryCell()) : - (DetId)(HGCHEDetId(id).geometryCell())); + DetId geoId; + if (topology().dddConstants().geomMode() == HGCalGeometryMode::Square) { + geoId = (id.subdetId() == HGCEE ? + (DetId)(HGCEEDetId(id).geometryCell()) : + (DetId)(HGCHEDetId(id).geometryCell())); + } else { + geoId = (DetId)(HGCalDetId(id).geometryCell()); + } cellIndex = topology().detId2denseGeomId(geoId); #ifdef DebugLog std::cout << "indexFor " << std::hex << id.rawId() << ":" << geoId.rawId() diff --git a/Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc b/Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc index 84950315a12ce..c8c86cfa0fd6e 100644 --- a/Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc +++ b/Geometry/HGCalGeometry/src/HGCalGeometryLoader.cc @@ -20,10 +20,11 @@ HGCalGeometry* HGCalGeometryLoader::build (const HGCalTopology& topology) { // allocate geometry HGCalGeometry* geom = new HGCalGeometry (topology); unsigned int numberOfCells = topology.totalGeomModules(); // both sides + unsigned int numberExpected= topology.allGeomModules(); #ifdef DebugLog - std::cout << "Number of Cells " << numberOfCells << " for sub-detector " - << topology.subDetector() << " Shape parameters " - << HGCalGeometry::k_NumberOfShapes << ":" + std::cout << "Number of Cells " << numberOfCells << ":" << numberExpected + << " for sub-detector " << topology.subDetector() + << " Shape parameters " << HGCalGeometry::k_NumberOfShapes << ":" << HGCalGeometry::k_NumberOfParametersPerShape << std::endl; #endif geom->allocateCorners( numberOfCells ) ; @@ -33,93 +34,137 @@ HGCalGeometry* HGCalGeometryLoader::build (const HGCalTopology& topology) { bool detType = topology.detectorType(); // loop over modules - std::vector::const_iterator trItr; - std::vector::const_iterator volItr; + std::vector::const_iterator trItr; + std::vector::const_iterator volItr; ParmVec params(HGCalGeometry::k_NumberOfParametersPerShape,0); unsigned int counter(0); for (trItr = topology.dddConstants().getFirstTrForm(); trItr != topology.dddConstants().getLastTrForm(); ++trItr) { int zside = trItr->zp; int layer = trItr->lay; - int sector = trItr->sec; - int subSec = (detType ? trItr->subsec : 0); - const HepGeom::Transform3D ht3d (trItr->hr, trItr->h3v); - DetId detId= ((subdet == HGCEE) ? - (DetId)(HGCEEDetId(subdet,zside,layer,sector,subSec,0)) : - (DetId)(HGCHEDetId(subdet,zside,layer,sector,subSec,0))); #ifdef DebugLog std::cout << "HGCalGeometryLoader:: Z:Layer:Sector:Subsector " << zside - << ":" << layer << ":" << sector << ":" << subSec << " transf " - << ht3d.getTranslation() << " and " << ht3d.getRotation(); + << ":" << layer <lay == layer) { - double alpha = ((detType && subSec == 0) ? -fabs(volItr->alpha) : - fabs(volItr->alpha)); - params[0] = volItr->dz; - params[1] = params[2] = 0; - params[3] = params[7] = volItr->h; - params[4] = params[8] = volItr->bl; - params[5] = params[9] = volItr->tl; - params[6] = params[10]= alpha; - params[11]= volItr->cellSize; + if (topology.geomMode() == HGCalGeometryMode::Square) { + int sector = trItr->sec; + int subSec = (detType ? trItr->subsec : 0); + const HepGeom::Transform3D ht3d (trItr->hr, trItr->h3v); + DetId detId= ((subdet == HGCEE) ? + (DetId)(HGCEEDetId(subdet,zside,layer,sector,subSec,0)) : + (DetId)(HGCHEDetId(subdet,zside,layer,sector,subSec,0))); #ifdef DebugLog - std::cout << "Volume Parameters"; - for (unsigned int i=0; i<12; ++i) std::cout << " : " << params[i]; - std::cout << std::endl; + std::cout << "HGCalGeometryLoader:: Sector:Subsector " << sector << ":" + << subSec << " transf " << ht3d.getTranslation() << " and " + << ht3d.getRotation(); #endif - std::vector corners (8); - - FlatTrd::createCorners( params, ht3d, corners ) ; - - const CCGFloat* parmPtr (CaloCellGeometry::getParmPtr(params, - geom->parMgr(), - geom->parVecVec() ) ) ; - - GlobalPoint front ( 0.25*( corners[0].x() + - corners[1].x() + - corners[2].x() + - corners[3].x() ), - 0.25*( corners[0].y() + - corners[1].y() + - corners[2].y() + - corners[3].y() ), - 0.25*( corners[0].z() + - corners[1].z() + - corners[2].z() + - corners[3].z() ) ) ; - - GlobalPoint back ( 0.25*( corners[4].x() + - corners[5].x() + - corners[6].x() + - corners[7].x() ), - 0.25*( corners[4].y() + - corners[5].y() + - corners[6].y() + - corners[7].y() ), - 0.25*( corners[4].z() + - corners[5].z() + - corners[6].z() + - corners[7].z() ) ) ; - - if (front.mag2() > back.mag2()) { // front should always point to the center, so swap front and back - std::swap (front, back); - std::swap_ranges (corners.begin(), corners.begin()+4, corners.begin()+4); + for (volItr = topology.dddConstants().getFirstModule(true); + volItr != topology.dddConstants().getLastModule(true); ++volItr) { + if (volItr->lay == layer) { + double alpha = ((detType && subSec == 0) ? -fabs(volItr->alpha) : + fabs(volItr->alpha)); + params[0] = volItr->dz; + params[1] = params[2] = 0; + params[3] = params[7] = volItr->h; + params[4] = params[8] = volItr->bl; + params[5] = params[9] = volItr->tl; + params[6] = params[10]= alpha; + params[11]= volItr->cellSize; + + buildGeom(params, ht3d, detId, geom); + counter++; + break; + } + } + } else { + for (int wafer=0; wafer 0) ? w.x() : -w.x(); + CLHEP::Hep3Vector h3v(xx,w.y(),trItr->h3v.z()); + const HepGeom::Transform3D ht3d (trItr->hr, h3v); +#ifdef DebugLog + std::cout << "HGCalGeometryLoader:: Wafer:Type " << wafer << ":" + << type << " transf " << ht3d.getTranslation() << " and " + << ht3d.getRotation(); +#endif + volItr = topology.dddConstants().getModule(wafer); + params[0] = volItr->dz; + params[1] = params[2] = 0; + params[3] = params[7] = volItr->h; + params[4] = params[8] = volItr->bl; + params[5] = params[9] = volItr->tl; + params[6] = params[10]= 0; + params[11]= topology.dddConstants().cellSizeHex(type); + + buildGeom(params, ht3d, detId, geom); + counter++; } - - geom->newCell(front, back, corners[0], parmPtr, detId) ; - counter++; - break; } } } + geom->sortDetIds(); - if (counter != numberOfCells) { - std::cerr << "inconsistent # of cells: expected " << numberOfCells << " , inited " << counter << std::endl; + if (counter != numberExpected) { + std::cerr << "inconsistent # of cells: expected " << numberExpected << ":" + << numberOfCells << " , inited " << counter << std::endl; assert( counter == numberOfCells ) ; } return geom; } + +void HGCalGeometryLoader::buildGeom(const ParmVec& params, + const HepGeom::Transform3D& ht3d, + const DetId& detId, HGCalGeometry* geom) { + +#ifdef DebugLog + std::cout << "Volume Parameters"; + for (unsigned int i=0; i<12; ++i) std::cout << " : " << params[i]; + std::cout << std::endl; +#endif + std::vector corners (8); + + FlatTrd::createCorners( params, ht3d, corners ) ; + + const CCGFloat* parmPtr (CaloCellGeometry::getParmPtr(params, + geom->parMgr(), + geom->parVecVec() ) ) ; + + GlobalPoint front ( 0.25*( corners[0].x() + + corners[1].x() + + corners[2].x() + + corners[3].x() ), + 0.25*( corners[0].y() + + corners[1].y() + + corners[2].y() + + corners[3].y() ), + 0.25*( corners[0].z() + + corners[1].z() + + corners[2].z() + + corners[3].z() ) ) ; + + GlobalPoint back ( 0.25*( corners[4].x() + + corners[5].x() + + corners[6].x() + + corners[7].x() ), + 0.25*( corners[4].y() + + corners[5].y() + + corners[6].y() + + corners[7].y() ), + 0.25*( corners[4].z() + + corners[5].z() + + corners[6].z() + + corners[7].z() ) ) ; + + if (front.mag2() > back.mag2()) { // front should always point to the center, so swap front and back + std::swap (front, back); + std::swap_ranges (corners.begin(), corners.begin()+4, corners.begin()+4); + } + + geom->newCell(front, back, corners[0], parmPtr, detId) ; +} diff --git a/Geometry/HGCalGeometry/test/HGCalGeometryTester.cc b/Geometry/HGCalGeometry/test/HGCalGeometryTester.cc index a76a8f0e8391f..56f85ba32e4c9 100644 --- a/Geometry/HGCalGeometry/test/HGCalGeometryTester.cc +++ b/Geometry/HGCalGeometry/test/HGCalGeometryTester.cc @@ -30,9 +30,15 @@ class HGCalGeometryTester : public edm::one::EDAnalyzer<> { private: void doTest(const HGCalGeometry& geom, ForwardSubdetector subdet); + + std::string name; + bool squareCell; }; -HGCalGeometryTester::HGCalGeometryTester(const edm::ParameterSet& ) {} +HGCalGeometryTester::HGCalGeometryTester(const edm::ParameterSet& iC) { + name = iC.getParameter("Detector"); + squareCell = iC.getParameter("SquareCell"); +} HGCalGeometryTester::~HGCalGeometryTester() {} @@ -40,27 +46,17 @@ HGCalGeometryTester::~HGCalGeometryTester() {} void HGCalGeometryTester::analyze(const edm::Event& , const edm::EventSetup& iSetup ) { - std::string name; - edm::ESHandle geom; - - name = "HGCalEESensitive"; - iSetup.get().get(name,geom); - if (geom.isValid()) doTest(*geom, HGCEE); - else std::cout << "Cannot get valid HGCalGeometry Object for " - << name << std::endl; + ForwardSubdetector subdet; + if (name == "HGCalHESiliconSensitive") subdet = HGCHEF; + else if (name == "HGCalHEScintillatorSensitive") subdet = HGCHEB; + else subdet = HGCEE; - name = "HGCalHESiliconSensitive"; + edm::ESHandle geom; iSetup.get().get(name,geom); - if (geom.isValid()) doTest(*geom, HGCHEF); - else std::cout << "Cannot get valid HGCalGeometry Object for " - << name << std::endl; - name = "HGCalHEScintillatorSensitive"; - iSetup.get().get(name,geom); - if (geom.isValid()) doTest(*geom,HGCHEB); + if (geom.isValid()) doTest(*geom, subdet); else std::cout << "Cannot get valid HGCalGeometry Object for " << name << std::endl; - } void HGCalGeometryTester::doTest(const HGCalGeometry& geom, @@ -73,18 +69,27 @@ void HGCalGeometryTester::doTest(const HGCalGeometry& geom, int sectors[]= {1, 7, 13}; int layers[] = {1, 5, 10}; int zsides[] = {1, -1}; - int cells[] = {1, 101, 201}; + int cells[] = {1, 51, 101}; + int wafers[] = {1, 101, 201, 301, 401}; + int types[] = {1, 0, 1, 0, 1}; + int ismax = (squareCell) ? 3 : 5; for (int iz = 0; iz < 2; ++iz) { int zside = zsides[iz]; - for (int is = 0; is < 3; ++is) { - int sector = sectors[is]; + for (int is = 0; is < ismax; ++is) { + int sector = (squareCell) ? sectors[is] : wafers[is]; + int type = (squareCell) ? 0 : types[is]; for (int il = 0; il < 3; ++il) { int layer = layers[il]; for (int ic = 0; ic < 3; ++ic) { int cell = cells[ic]; - DetId id1= ((subdet == HGCEE) ? - (DetId)(HGCEEDetId(subdet,zside,layer,sector,0,cell)) : - (DetId)(HGCHEDetId(subdet,zside,layer,sector,0,cell))); + DetId id1; + if (squareCell) { + id1 = ((subdet == HGCEE) ? + (DetId)(HGCEEDetId(subdet,zside,layer,sector,type,cell)) : + (DetId)(HGCHEDetId(subdet,zside,layer,sector,type,cell))); + } else { + id1 = (DetId)(HGCalDetId(subdet,zside,layer,type,sector,cell)); + } const CaloCellGeometry* icell1 = geom.getGeometry(id1); GlobalPoint global1 = geom.getPosition(id1); DetId idc1 = geom.getClosestCell(global1); @@ -92,39 +97,51 @@ void HGCalGeometryTester::doTest(const HGCalGeometry& geom, << ":" << sector << ":0:" << cell << ") Geom " << icell1 << " position (" << global1.x() << ", " << global1.y() << ", " << global1.z() << ") ids " << std::hex - << id1.rawId() << ":" << idc1.rawId() << std::dec - << " parameter[11] = " << icell1->param()[10] << ":" + << id1.rawId() << ":" << idc1.rawId() << std::dec; + if (squareCell) { + if (subdet == HGCEE) + std::cout << ":" << HGCEEDetId(id1) << ":" << HGCEEDetId(idc1); + else + std::cout << ":" << HGCHEDetId(id1) << ":" << HGCHEDetId(idc1); + } else { + std::cout << ":" << HGCalDetId(id1) << ":" << HGCalDetId(idc1); + } + std::cout << " parameter[11] = " << icell1->param()[10] << ":" << icell1->param()[11] << std::endl; if (id1.rawId() != idc1.rawId()) std::cout << "***** ERROR *****\n"; - DetId id2= ((subdet == HGCEE) ? - (DetId)(HGCEEDetId(subdet,zside,layer,sector,1,cell)) : - (DetId)(HGCHEDetId(subdet,zside,layer,sector,1,cell))); - const CaloCellGeometry* icell2 = geom.getGeometry(id2); - GlobalPoint global2 = geom.getPosition(id2); - DetId idc2 = geom.getClosestCell(global2); - std::cout << "DetId (" << subdet << ":" << zside << ":" << layer - << ":" << sector << ":1:" << cell << ") Geom " << icell2 - << " position (" << global2.x() << ", " << global2.y() - << ", " << global2.z() << ") ids " << std::hex - << id2.rawId() << ":" << idc2.rawId() << std::dec - << " parameter[11] = " << icell2->param()[10] << ":" - << icell2->param()[11] << std::endl; - if (id2.rawId() != idc2.rawId()) std::cout << "***** ERROR *****\n"; + if (squareCell) { + DetId id2= ((subdet == HGCEE) ? + (DetId)(HGCEEDetId(subdet,zside,layer,sector,1,cell)) : + (DetId)(HGCHEDetId(subdet,zside,layer,sector,1,cell))); + const CaloCellGeometry* icell2 = geom.getGeometry(id2); + GlobalPoint global2 = geom.getPosition(id2); + DetId idc2 = geom.getClosestCell(global2); + std::cout << "DetId (" << subdet << ":" << zside << ":" << layer + << ":" << sector << ":1:" << cell << ") Geom " << icell2 + << " position (" << global2.x() << ", " << global2.y() + << ", " << global2.z() << ") ids " << std::hex + << id2.rawId() << ":" << idc2.rawId() << std::dec + << " parameter[11] = " << icell2->param()[10] << ":" + << icell2->param()[11] << std::endl; + if (id2.rawId() != idc2.rawId()) std::cout << "***** ERROR *****\n"; + } } } } } - uint32_t probids[] = {1711603886, 1711603890, 1761408735, 1761411303, - 1801744385, 1805447194}; - for (int k=0; k<6; ++k) { - DetId id(probids[k]); - if (id.det() == DetId::Forward && id.subdetId() == (int)(subdet)) { - if (subdet == HGCEE) std::cout << "Test " << HGCEEDetId(id) << std::endl; - else std::cout << "Test " << HGCHEDetId(id) << std::endl; - const CaloCellGeometry* icell = geom.getGeometry(id); - GlobalPoint global = geom.getPosition(id); - std::cout << "Geom Cell: " << icell << " position (" << global.x() - << ", " << global.y() << ", " << global.z() << ")"<< std::endl; + if (squareCell) { + uint32_t probids[] = {1711603886, 1711603890, 1761408735, 1761411303, + 1801744385, 1805447194}; + for (int k=0; k<6; ++k) { + DetId id(probids[k]); + if (id.det() == DetId::Forward && id.subdetId() == (int)(subdet)) { + if (subdet == HGCEE) std::cout << "Test " << HGCEEDetId(id) << std::endl; + else std::cout << "Test " << HGCHEDetId(id) << std::endl; + const CaloCellGeometry* icell = geom.getGeometry(id); + GlobalPoint global = geom.getPosition(id); + std::cout << "Geom Cell: " << icell << " position (" << global.x() + << ", " << global.y() << ", " << global.z() << ")"<< std::endl; + } } } } diff --git a/Geometry/HGCalGeometry/test/python/testHGCalGeometryFromLocalDB_cfg.py b/Geometry/HGCalGeometry/test/python/testHGCalGeometryFromLocalDB_cfg.py index 7ca7fb8906bd4..496472b538081 100644 --- a/Geometry/HGCalGeometry/test/python/testHGCalGeometryFromLocalDB_cfg.py +++ b/Geometry/HGCalGeometry/test/python/testHGCalGeometryFromLocalDB_cfg.py @@ -39,6 +39,19 @@ input = cms.untracked.int32(1) ) -process.test = cms.EDAnalyzer("HGCalGeometryTester") +process.prodEE = cms.EDAnalyzer("HGCalGeometryTester", + Detector = cms.string("HGCalEESensitive"), + SquareCell = cms.bool(True), + ) + +process.prodHEF = process.prodEE.clone( + Detector = "HGCalHESiliconSensitive", + SquareCell = True +) + +process.prodHEB = process.prodEE.clone( + Detector = "HGCalHEScintillatorSensitive", + SquareCell = True +) -process.p1 = cms.Path(process.test) +process.p1 = cms.Path(process.prodEE+process.prodHEF+process.prodHEB) diff --git a/Geometry/HGCalGeometry/test/python/testHGCalGeometry_cfg.py b/Geometry/HGCalGeometry/test/python/testHGCalGeometry_cfg.py index d2230f5fa8e5b..f580e810995e8 100644 --- a/Geometry/HGCalGeometry/test/python/testHGCalGeometry_cfg.py +++ b/Geometry/HGCalGeometry/test/python/testHGCalGeometry_cfg.py @@ -3,8 +3,9 @@ process = cms.Process("PROD") process.load("SimGeneral.HepPDTESSource.pdt_cfi") -process.load("Geometry.HGCalCommonData.testHGCalXML_cfi") +process.load("Geometry.HGCalCommonData.testHGCXML_cfi") process.load("Geometry.HGCalCommonData.hgcalNumberingInitialization_cfi") +process.load("Geometry.HGCalCommonData.hgcalParametersInitialization_cfi") process.load("Geometry.CaloEventSetup.HGCalTopology_cfi") process.load("Geometry.HGCalGeometry.HGCalGeometryESProducer_cfi") @@ -47,6 +48,19 @@ input = cms.untracked.int32(1) ) -process.prod = cms.EDAnalyzer("HGCalGeometryTester") +process.prodEE = cms.EDAnalyzer("HGCalGeometryTester", + Detector = cms.string("HGCalEESensitive"), + SquareCell = cms.bool(False), + ) -process.p1 = cms.Path(process.generator*process.prod) +process.prodHEF = process.prodEE.clone( + Detector = "HGCalHESiliconSensitive", + SquareCell = False +) + +process.prodHEB = process.prodEE.clone( + Detector = "HGCalHEScintillatorSensitive", + SquareCell = True +) + +process.p1 = cms.Path(process.generator*process.prodEE*process.prodHEF) diff --git a/SimCalorimetry/HGCalSimProducers/src/HGCDigitizer.cc b/SimCalorimetry/HGCalSimProducers/src/HGCDigitizer.cc index 7f013b23c2a6f..2cfa3964b6365 100644 --- a/SimCalorimetry/HGCalSimProducers/src/HGCDigitizer.cc +++ b/SimCalorimetry/HGCalSimProducers/src/HGCDigitizer.cc @@ -13,6 +13,7 @@ #include "FWCore/Framework/interface/Event.h" #include "FWCore/Framework/interface/EventSetup.h" #include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h" #include #include @@ -170,76 +171,86 @@ void HGCDigitizer::accumulate(edm::Handle const &hits, //create list of tuples (pos in container, RECO DetId, time) to be sorted first int nchits=(int)hits->size(); std::vector< HGCCaloHitTuple_t > hitRefs(nchits); - for(int i=0; iat(i).id(); + for(int i=0; iat(i).id(); + if (dddConst.geomMode() == HGCalGeometryMode::Square) { HGCalTestNumbering::unpackSquareIndex(simId, zp, layer, sec, subsec, cell); - //skip this hit if after ganging it is not valid - std::pair recoLayerCell=dddConst.simToReco(cell,layer,topo.detectorType()); - cell = recoLayerCell.first; - layer = recoLayerCell.second; - if (layer<0 || cell<0) { - hitRefs[i]=std::make_tuple( i, 0, 0.); - continue; - } + } else { + int subdet; + HGCalTestNumbering::unpackHexagonIndex(simId, subdet, zp, layer, sec, subsec, cell); + mySubDet_ = (ForwardSubdetector)(subdet); + //sec is wafer and subsec is celltyp + } + //skip this hit if after ganging it is not valid + std::pair recoLayerCell=dddConst.simToReco(cell,layer,sec,topo.detectorType()); + cell = recoLayerCell.first; + layer = recoLayerCell.second; + if (layer<0 || cell<0) { + hitRefs[i]=std::make_tuple( i, 0, 0.); + continue; + } - //assign the RECO DetId - DetId id( producesEEDigis() ? - (uint32_t)HGCEEDetId(mySubDet_,zp,layer,sec,subsec,cell): - (uint32_t)HGCHEDetId(mySubDet_,zp,layer,sec,subsec,cell) ); + //assign the RECO DetId + DetId id; + if (dddConst.geomMode() == HGCalGeometryMode::Square) { + id = (producesEEDigis() ? + (uint32_t)HGCEEDetId(mySubDet_,zp,layer,sec,subsec,cell): + (uint32_t)HGCHEDetId(mySubDet_,zp,layer,sec,subsec,cell) ); + } else { + id = HGCalDetId(mySubDet_,zp,layer,subsec,sec,cell); + } - if (verbosity_>0) { - if (producesEEDigis()) + if (verbosity_>0) { + if (producesEEDigis()) edm::LogInfo("HGCDigitizer") <<" i/p " << std::hex << simId << std::dec << " o/p " << HGCEEDetId(id) << std::endl; - else + else edm::LogInfo("HGCDigitizer") << " i/p " << std::hex << simId << std::dec << " o/p " << HGCHEDetId(id) << std::endl; - } - - hitRefs[i]=std::make_tuple( i, - id.rawId(), - (float)hits->at(i).time() ); } + + hitRefs[i]=std::make_tuple( i, + id.rawId(), + (float)hits->at(i).time() ); + } std::sort(hitRefs.begin(),hitRefs.end(),this->orderByDetIdThenTime); //loop over sorted hits - for(int i=0; i(hitRefs[i]); - const uint32_t id = std::get<1>(hitRefs[i]); - if(id==0) continue; // to be ignored at RECO level + for(int i=0; i(hitRefs[i]); + const uint32_t id = std::get<1>(hitRefs[i]); + if(id==0) continue; // to be ignored at RECO level - const float toa = std::get<2>(hitRefs[i]); - const PCaloHit &hit=hits->at( hitidx ); - const float charge = hit.energy()*1e6*keV2fC; + const float toa = std::get<2>(hitRefs[i]); + const PCaloHit &hit=hits->at( hitidx ); + const float charge = hit.energy()*1e6*keV2fC; - //distance to the center of the detector - const float dist2center( geom->getPosition(id).mag() ); + //distance to the center of the detector + const float dist2center( geom->getPosition(id).mag() ); - //hit time: [time()]=ns [centerDist]=cm [refSpeed_]=cm/ns + delay by 1ns - //accumulate in 15 buckets of 25ns (9 pre-samples, 1 in-time, 5 post-samples) - const float tof = toa-dist2center/refSpeed_+tofDelay_ ; - const int itime= std::floor( tof/bxTime_ ) + 9; + //hit time: [time()]=ns [centerDist]=cm [refSpeed_]=cm/ns + delay by 1ns + //accumulate in 15 buckets of 25ns (9 pre-samples, 1 in-time, 5 post-samples) + const float tof = toa-dist2center/refSpeed_+tofDelay_ ; + const int itime= std::floor( tof/bxTime_ ) + 9; - //no need to add bx crossing - tof comes already corrected from the mixing module - //itime += bxCrossing; - //itime += 9; + //no need to add bx crossing - tof comes already corrected from the mixing module + //itime += bxCrossing; + //itime += 9; - if(itime<0 || itime>14) continue; - - //check if already existing (perhaps could remove this in the future - 2nd event should have all defined) - HGCSimHitDataAccumulator::iterator simHitIt=simHitAccumulator_->find(id); - if(simHitIt == simHitAccumulator_->end()) { - simHitIt = simHitAccumulator_->insert( std::make_pair(id,baseData) ).first; - } + if(itime<0 || itime>14) continue; + + //check if already existing (perhaps could remove this in the future - 2nd event should have all defined) + HGCSimHitDataAccumulator::iterator simHitIt=simHitAccumulator_->find(id); + if(simHitIt == simHitAccumulator_->end()) { + simHitIt = simHitAccumulator_->insert( std::make_pair(id,baseData) ).first; + } - //check if time index is ok and store energy - if(itime >= (int)simHitIt->second[0].size() ) continue; + //check if time index is ok and store energy + if(itime >= (int)simHitIt->second[0].size() ) continue; - (simHitIt->second)[0][itime] += charge; - float accCharge=(simHitIt->second)[0][itime]; + (simHitIt->second)[0][itime] += charge; + float accCharge=(simHitIt->second)[0][itime]; - //time-of-arrival (check how to be used) + //time-of-arrival (check how to be used) if(weightToAbyEnergy) (simHitIt->second)[1][itime] += charge*tof; else if((simHitIt->second)[1][itime]==0) { diff --git a/SimDataFormats/CaloTest/src/HGCalTestNumbering.cc b/SimDataFormats/CaloTest/src/HGCalTestNumbering.cc index 8cae7c3bb83b4..69c245f809424 100644 --- a/SimDataFormats/CaloTest/src/HGCalTestNumbering.cc +++ b/SimDataFormats/CaloTest/src/HGCalTestNumbering.cc @@ -1,22 +1,12 @@ #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h" #include -//#define DebugLog +#define DebugLog uint32_t HGCalTestNumbering::packSquareIndex(int zp, int lay, int sec, int subsec, int cell) { - if (cell > kHGCalCellSMask || sec>kHGCalSectorSMask || - subsec > kHGCalSubSectorSMask || lay>kHGCalLayerSMask ) { -#ifdef DebugLog - std::cout << "[HGCalTestNumbering] request for new id for layer=" << lay - << " zp=" << zp - << " sector=" << sec - << " subsec=" << subsec - << " cell=" << cell - << " has one or more fields out of bounds and will be reset" - << std::endl; -#endif + if (!HGCalTestNumbering::isValidSquare(zp, lay, sec, subsec,lay)) { zp = lay = sec = subsec = cell = 0; } @@ -34,19 +24,7 @@ uint32_t HGCalTestNumbering::packHexagonIndex(int subdet, int zp, int lay, int wafer, int celltyp, int cell) { - if (cell > kHGCalCellHMask || celltyp>kHGCalCellTypHMask || - wafer > kHGCalWaferHMask || lay>kHGCalLayerSMask || - subdet > kHGCalSubdetHMask) { -#ifdef DebugLog - std::cout << "[HGCalTestNumbering] request for new id for layer=" << lay - << " zp=" << zp - << " wafer=" << wafer - << " celltyp=" << celltyp - << " cell=" << cell - << " for subdet=" << subdet - << " has one or more fields out of bounds and will be reset" - << std::endl; -#endif + if (!HGCalTestNumbering::isValidHexagon(subdet,zp,lay,wafer,celltyp,cell)) { subdet = zp = lay = wafer = celltyp = cell = 0; } diff --git a/SimG4CMS/Calo/interface/HGCNumberingScheme.h b/SimG4CMS/Calo/interface/HGCNumberingScheme.h index 8f953818a6bc5..9673a2390a36e 100644 --- a/SimG4CMS/Calo/interface/HGCNumberingScheme.h +++ b/SimG4CMS/Calo/interface/HGCNumberingScheme.h @@ -21,15 +21,16 @@ class HGCNumberingScheme : public CaloNumberingScheme { enum HGCNumberingParameters { HGCCellSize }; - HGCNumberingScheme(const DDCompactView & cpv, std::string& name, bool check, - int verbose); + HGCNumberingScheme(HGCalDDDConstants* hgc, std::string& name, + bool check, int verbose); virtual ~HGCNumberingScheme(); /** @short assigns the det id to a hit */ - virtual uint32_t getUnitID(ForwardSubdetector subdet, int layer, int module, int iz, const G4ThreeVector &pos); + uint32_t getUnitID(ForwardSubdetector subdet, int layer, int module, + int cell, int iz, const G4ThreeVector &pos); /** @short maps a hit position to a sequential cell in a trapezoid surface defined by h,b,t diff --git a/SimG4CMS/Calo/interface/HGCSD.h b/SimG4CMS/Calo/interface/HGCSD.h index edb433f948681..04c95ee555666 100644 --- a/SimG4CMS/Calo/interface/HGCSD.h +++ b/SimG4CMS/Calo/interface/HGCSD.h @@ -7,6 +7,8 @@ /////////////////////////////////////////////////////////////////////////////// #include "SimG4CMS/Calo/interface/CaloSD.h" +#include "SimG4Core/Notification/interface/BeginOfJob.h" +#include "Geometry/HGCalCommonData/interface/HGCalDDDConstants.h" #include "SimG4CMS/Calo/interface/HGCNumberingScheme.h" #include "DetectorDescription/Core/interface/DDsvalues.h" @@ -21,7 +23,7 @@ class G4LogicalVolume; class G4Material; class G4Step; -class HGCSD : public CaloSD { +class HGCSD : public CaloSD, public Observer { public: @@ -34,15 +36,19 @@ class HGCSD : public CaloSD { protected: + virtual void update(const BeginOfJob *); virtual void initRun(); virtual bool filterHit(CaloG4Hit*, double); private: - uint32_t setDetUnitId(ForwardSubdetector &, int &, int &, int &, G4ThreeVector &); + uint32_t setDetUnitId(ForwardSubdetector&, int, int, int, int, G4ThreeVector &); bool isItinFidVolume (G4ThreeVector&) {return true;} int setTrackID(G4Step * step); + std::string nameX; + bool checkID; + HGCalDDDConstants* hgcalCons; HGCNumberingScheme* numberingScheme; G4int verbosity, mumPDG, mupPDG; double eminHit; diff --git a/SimG4CMS/Calo/src/HGCNumberingScheme.cc b/SimG4CMS/Calo/src/HGCNumberingScheme.cc index 070762e772525..411735d10bca4 100644 --- a/SimG4CMS/Calo/src/HGCNumberingScheme.cc +++ b/SimG4CMS/Calo/src/HGCNumberingScheme.cc @@ -4,6 +4,7 @@ /////////////////////////////////////////////////////////////////////////////// #include "SimG4CMS/Calo/interface/HGCNumberingScheme.h" #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h" +#include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h" #include "DataFormats/Math/interface/FastMath.h" #include "CLHEP/Units/GlobalSystemOfUnits.h" @@ -11,11 +12,11 @@ //#define DebugLog -HGCNumberingScheme::HGCNumberingScheme(const DDCompactView & cpv, +HGCNumberingScheme::HGCNumberingScheme(HGCalDDDConstants* hgc, std::string & name, bool check, int verbose) : CaloNumberingScheme(0), check_(check), verbosity(verbose), - hgcons(new HGCalDDDConstants(cpv,name)) { + hgcons(hgc) { edm::LogInfo("HGCSim") << "Creating HGCNumberingScheme for " << name; } @@ -24,34 +25,60 @@ HGCNumberingScheme::~HGCNumberingScheme() { } // -uint32_t HGCNumberingScheme::getUnitID(ForwardSubdetector subdet, int layer, int sector, int iz, const G4ThreeVector &pos) { +uint32_t HGCNumberingScheme::getUnitID(ForwardSubdetector subdet, int layer, + int module, int cell, int iz, + const G4ThreeVector &pos) { + // module is the sector # for square cell + // the copy number of the wafer as placed in the layer + int phiSector(0), icell(0), celltyp(0), wafer(0); + uint32_t index(0); + if (hgcons->geomMode() == HGCalGeometryMode::Square) { + std::pair phicell = hgcons->assignCell(pos.x(),pos.y(),layer,0,false); + phiSector = phicell.first; + icell = phicell.second; - std::pair phicell = hgcons->assignCell(pos.x(),pos.y(),layer,0,false); - int phiSector = phicell.first; - int icell = phicell.second; - - //build the index - uint32_t index = HGCalTestNumbering::packSquareIndex(iz,layer,sector,phiSector,icell); - - //check if it fits - if ((!HGCalTestNumbering::isValidSquare(iz,layer,sector,phiSector,icell)) || - (!hgcons->isValid(layer,sector,icell,false))) { - index = 0; - if (check_ && icell != -1) { - edm::LogError("HGCSim") << "[HGCNumberingScheme] ID out of bounds :" - << " Subdet= " << subdet << " Zside= " << iz - << " Layer= " << layer << " Sector= " << sector - << " SubSector= " << phiSector << " Cell= " - << icell << " Local position: (" << pos.x() - << "," << pos.y() << "," << pos.z() << ")"; + //build the index + index = HGCalTestNumbering::packSquareIndex(iz,layer,module,phiSector,icell); + //check if it fits + if (!hgcons->isValid(layer,module,icell,false)) { + index = 0; + if (check_ && icell != -1) + edm::LogError("HGCSim") << "[HGCNumberingScheme] ID out of bounds :" + << " Subdet= " << subdet << " Zside= " << iz + << " Layer= " << layer << " Module= " << module + << " SubSector= " << phiSector << " Cell= " + << icell << " Local position: (" << pos.x() + << "," << pos.y() << "," << pos.z() << ")"; + } + } else { + celltyp = cell/1000; + icell = cell%1000; + if (celltyp != 1) celltyp = 0; + wafer = hgcons->waferFromCopy(module); + index = HGCalTestNumbering::packHexagonIndex((int)subdet,iz,layer,wafer, + celltyp,icell); + //check if it fits + if (!hgcons->isValid(layer,wafer,icell,false)) { + index = 0; + if (check_) { + edm::LogError("HGCSim") << "[HGCNumberingScheme] ID out of bounds :" + << " Subdet= " << subdet << " Zside= " << iz + << " Layer= " << layer << " Wafer= " << module + << " CellType= " << celltyp << " Cell= " + << icell; + } } - } + } #ifdef DebugLog - if (verbosity > 0) + if (verbosity > 0) { std::cout << "HGCNumberingScheme::i/p " << subdet << ":" << layer << ":" - << sector << ":" << iz << ":" << pos << " o/p " << phiSector - << ":" << icell << ":" << std::hex << index << std::dec - << std::endl; + << module << ":" << iz << ":"; + if (hgcons->geomMode() == HGCalGeometryMode::Square) + std::cout << pos << " o/p " << phiSector << ":" << icell; + else + std::cout << wafer << ":" << celltyp << ":" << icell; + std::cout << ":" << std::hex << index << std::dec << std::endl; + } #endif return index; } diff --git a/SimG4CMS/Calo/src/HGCSD.cc b/SimG4CMS/Calo/src/HGCSD.cc index 0b847d41358b1..afd72e12e0707 100644 --- a/SimG4CMS/Calo/src/HGCSD.cc +++ b/SimG4CMS/Calo/src/HGCSD.cc @@ -13,7 +13,10 @@ #include "DetectorDescription/Core/interface/DDMaterial.h" #include "DetectorDescription/Core/interface/DDValue.h" #include "FWCore/Utilities/interface/Exception.h" - +#include "FWCore/Framework/interface/ESHandle.h" +#include "FWCore/Framework/interface/EventSetup.h" +#include "Geometry/Records/interface/IdealGeometryRecord.h" +#include "Geometry/HGCalCommonData/interface/HGCalGeometryMode.h" #include "G4LogicalVolumeStore.hh" #include "G4LogicalVolume.hh" #include "G4Step.hh" @@ -34,17 +37,17 @@ HGCSD::HGCSD(G4String name, const DDCompactView & cpv, CaloSD(name, cpv, clg, p, manager, (float)(p.getParameter("HGCSD").getParameter("TimeSliceUnit")), p.getParameter("HGCSD").getParameter("IgnoreTrackID")), - numberingScheme(0) { + hgcalCons(0), numberingScheme(0) { edm::ParameterSet m_HGC = p.getParameter("HGCSD"); eminHit = m_HGC.getParameter("EminHit")*CLHEP::MeV; - bool checkID = m_HGC.getUntrackedParameter("CheckID", false); + checkID = m_HGC.getUntrackedParameter("CheckID", false); verbosity = m_HGC.getUntrackedParameter("Verbosity",0); //this is defined in the hgcsens.xml G4String myName(this->nameOfSD()); myFwdSubdet_= ForwardSubdetector::ForwardEmpty; - std::string nameX("HGCal"); + nameX = "HGCal"; if (myName.find("HitsEE")!=std::string::npos) { myFwdSubdet_ = ForwardSubdetector::HGCEE; nameX = "HGCalEESensitive"; @@ -67,8 +70,6 @@ HGCSD::HGCSD(G4String name, const DDCompactView & cpv, << "**************************************************"; #endif edm::LogInfo("HGCSim") << "HGCSD:: Threshold for storing hits: " << eminHit; - - numberingScheme = new HGCNumberingScheme(cpv,nameX,checkID,verbosity); } HGCSD::~HGCSD() { @@ -121,12 +122,36 @@ uint32_t HGCSD::setDetUnitId(G4Step * aStep) { //get the det unit id with ForwardSubdetector subdet = myFwdSubdet_; - int layer = touch->GetReplicaNumber(0); - int module = touch->GetReplicaNumber(1); + int layer(0), module(0), cell(0); + if (hgcalCons->geomMode() == HGCalGeometryMode::Square) { + layer = touch->GetReplicaNumber(0); + module = touch->GetReplicaNumber(1); + } else { + layer = touch->GetReplicaNumber(2); + module = touch->GetReplicaNumber(1); + cell = touch->GetReplicaNumber(0); + } if (verbosity > 0) - std::cout << "HGCSD::Global " << hitPoint << " local " << localpos - << std::endl; - return setDetUnitId (subdet, layer, module, iz, localpos); + edm::LogInfo("HGCSim") << "HGCSD::Global " << hitPoint << " local " + << localpos << "Layer|Module|Cell " << layer << "|" + << module << "|" << cell; + return setDetUnitId (subdet, layer, module, cell, iz, localpos); +} + +void HGCSD::update(const BeginOfJob * job) { + + const edm::EventSetup* es = (*job)(); + edm::ESHandle hdc; + es->get().get(nameX,hdc); + if (hdc.isValid()) { + hgcalCons = (HGCalDDDConstants*)(&(*hdc)); + } else { + edm::LogError("HGCSim") << "HCalSD : Cannot find HGCalDDDConstants for " + << nameX; + throw cms::Exception("Unknown", "HGCSD") << "Cannot find HGCalDDDConstants for " << nameX << "\n"; + } + + numberingScheme = new HGCNumberingScheme(hgcalCons,nameX,checkID,verbosity); } void HGCSD::initRun() { @@ -146,8 +171,9 @@ bool HGCSD::filterHit(CaloG4Hit* aHit, double time) { // -uint32_t HGCSD::setDetUnitId (ForwardSubdetector &subdet, int &layer, int &module, int &iz, G4ThreeVector &pos) { - return (numberingScheme ? numberingScheme->getUnitID(subdet, layer, module, iz, pos) : 0); +uint32_t HGCSD::setDetUnitId (ForwardSubdetector &subdet, int layer, int module, + int cell, int iz, G4ThreeVector &pos) { + return (numberingScheme ? numberingScheme->getUnitID(subdet, layer, module, cell, iz, pos) : 0); } // From c1209d45593b475a3844ad15d7fba2e67cdab314 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Sun, 27 Dec 2015 01:54:42 +0100 Subject: [PATCH 2/6] Take into account Iana's comments --- .../interface/HGCalParameters.h | 63 +++++- .../GeometryObjects/src/HGCalParameters.cc | 146 +++++++++++++ .../GeometryObjects/src/classes_def.xml | 46 ++++- Geometry/CaloGeometry/src/FlatTrd.cc | 10 + Geometry/HGCalCommonData/BuildFile.xml | 1 + .../interface/HGCalDDDConstants.h | 19 +- .../plugins/DDHGCalWaferAlgo.cc | 2 +- .../HGCalCommonData/src/HGCalDDDConstants.cc | 153 ++++++++------ .../src/HGCalGeomParameters.cc | 193 ++++++++++-------- .../test/HGCalNumberingTester.cc | 4 +- .../plugins/HGCalGeometryESProducer.cc | 2 +- Geometry/HGCalGeometry/src/HGCalGeometry.cc | 78 ++++--- .../HGCalGeometry/src/HGCalGeometryLoader.cc | 61 +++--- .../HGCalGeometry/test/HGCalGeometryTester.cc | 4 +- .../CaloTest/src/HGCalTestNumbering.cc | 2 +- 15 files changed, 536 insertions(+), 248 deletions(-) create mode 100644 CondFormats/GeometryObjects/src/HGCalParameters.cc diff --git a/CondFormats/GeometryObjects/interface/HGCalParameters.h b/CondFormats/GeometryObjects/interface/HGCalParameters.h index 5650937e412e8..1d7769cc03bc7 100644 --- a/CondFormats/GeometryObjects/interface/HGCalParameters.h +++ b/CondFormats/GeometryObjects/interface/HGCalParameters.h @@ -6,14 +6,12 @@ #include #include #include +#include #include class HGCalParameters { public: - - HGCalParameters(const std::string& nam): name_(nam) { } - ~HGCalParameters( void ) { } struct hgtrap { int lay; @@ -22,18 +20,59 @@ class HGCalParameters { struct hgtrform { int zp, lay, sec, subsec; - bool used; CLHEP::Hep3Vector h3v; CLHEP::HepRotation hr; }; + + HGCalParameters(const std::string& nam); + ~HGCalParameters( void ); + void fillModule(const hgtrap& mytr, bool reco); + hgtrap getModule(unsigned int k, bool reco) const; + void fillTrForm(const hgtrform& mytr); + hgtrform getTrForm(unsigned int k) const; + void addTrForm(const CLHEP::Hep3Vector& h3v); + void scaleTrForm(double); + + static const int kMaskZside = 0x1; + static const int kMaskLayer = 0x7F; + static const int kMaskSector = 0x3FF; + static const int kMaskSubSec = 0x1; + static const int kShiftZside = 19; + static const int kShiftLayer = 12; + static const int kShiftSector = 1; + static const int kShiftSubSec = 0; std::string name_; int nCells_; int nSectors_; std::vector cellSize_; - std::vector modules_; - std::vector moduler_; - std::vector trform_; + std::vector moduleLayS_; + std::vector moduleBlS_; + std::vector moduleTlS_; + std::vector moduleHS_; + std::vector moduleDzS_; + std::vector moduleAlphaS_; + std::vector moduleCellS_; + std::vector moduleLayR_; + std::vector moduleBlR_; + std::vector moduleTlR_; + std::vector moduleHR_; + std::vector moduleDzR_; + std::vector moduleAlphaR_; + std::vector moduleCellR_; + std::vector trformIndex_; + std::vector trformTranX_; + std::vector trformTranY_; + std::vector trformTranZ_; + std::vector trformRotXX_; + std::vector trformRotYX_; + std::vector trformRotZX_; + std::vector trformRotXY_; + std::vector trformRotYY_; + std::vector trformRotZY_; + std::vector trformRotXZ_; + std::vector trformRotYZ_; + std::vector trformRotZZ_; std::vector layer_; std::vector layerIndex_; std::vector layerGroup_; @@ -47,15 +86,19 @@ class HGCalParameters { std::vector waferCopy_; std::vector waferTypeL_; std::vector waferTypeT_; - std::vector waferPos_; - std::vector cellFine_; - std::vector cellCoarse_; + std::vector waferPosX_; + std::vector waferPosY_; + std::vector cellFineX_; + std::vector cellFineY_; + std::vector cellCoarseX_; + std::vector cellCoarseY_; std::vector layerGroupM_; std::vector layerGroupO_; std::vector boundR_; double waferR_; int mode_; + COND_SERIALIZABLE; }; #endif diff --git a/CondFormats/GeometryObjects/src/HGCalParameters.cc b/CondFormats/GeometryObjects/src/HGCalParameters.cc new file mode 100644 index 0000000000000..ab4c7c7729919 --- /dev/null +++ b/CondFormats/GeometryObjects/src/HGCalParameters.cc @@ -0,0 +1,146 @@ +#include "CondFormats/GeometryObjects/interface/HGCalParameters.h" +//#define DebugLog + +HGCalParameters::HGCalParameters(const std::string& nam): name_(nam) { } + +HGCalParameters::~HGCalParameters() { } + +void HGCalParameters::fillModule(const HGCalParameters::hgtrap& mytr, + bool reco) { + + if (reco) { + moduleLayR_.push_back(mytr.lay); + moduleBlR_.push_back(mytr.bl); + moduleTlR_.push_back(mytr.tl); + moduleHR_.push_back(mytr.h); + moduleDzR_.push_back(mytr.dz); + moduleAlphaR_.push_back(mytr.alpha); + moduleCellR_.push_back(mytr.cellSize); + } else { + moduleLayS_.push_back(mytr.lay); + moduleBlS_.push_back(mytr.bl); + moduleTlS_.push_back(mytr.tl); + moduleHS_.push_back(mytr.h); + moduleDzS_.push_back(mytr.dz); + moduleAlphaS_.push_back(mytr.alpha); + moduleCellS_.push_back(mytr.cellSize); + } +} + +HGCalParameters::hgtrap HGCalParameters::getModule(unsigned int k, + bool reco) const { + HGCalParameters::hgtrap mytr; + if (reco) { + if (k < moduleLayR_.size()) { + mytr.lay = moduleLayR_[k]; + mytr.bl = moduleBlR_[k]; + mytr.tl = moduleTlR_[k]; + mytr.h = moduleHR_[k]; + mytr.dz = moduleDzR_[k]; + mytr.alpha = moduleAlphaR_[k]; + mytr.cellSize = moduleCellR_[k]; + } else { + mytr.lay = -1; + mytr.bl = mytr.tl = mytr.h = mytr.dz = mytr.alpha = mytr.cellSize = 0; + } + } else { + if (k < moduleLayS_.size()) { + mytr.lay = moduleLayS_[k]; + mytr.bl = moduleBlS_[k]; + mytr.tl = moduleTlS_[k]; + mytr.h = moduleHS_[k]; + mytr.dz = moduleDzS_[k]; + mytr.alpha = moduleAlphaS_[k]; + mytr.cellSize = moduleCellS_[k]; + } else { + mytr.lay = -1; + mytr.bl = mytr.tl = mytr.h = mytr.dz = mytr.alpha = mytr.cellSize = 0; + } + } + return mytr; +} + +void HGCalParameters::fillTrForm(const HGCalParameters::hgtrform& mytr) { + + int zp = (mytr.zp == 1) ? 1 : 0; + uint32_t indx = ((zp & kMaskZside) << kShiftZside); + indx |= ((mytr.lay & kMaskLayer) << kShiftLayer); + indx |= ((mytr.sec & kMaskSector) << kShiftSector); + indx |= ((mytr.subsec & kMaskSubSec) << kShiftSubSec); +// std::cout << "ZP " << zp << ":" << kMaskZside << ":" << kShiftZside << ((zp & kMaskZside) << kShiftZside) << " Lay " << mytr.lay << ":" << kMaskLayer << ":" << kShiftLayer << ":" << ((mytr.lay & kMaskLayer) << kShiftLayer) << " Sector " << mytr.sec << ":" << kMaskSector << ":" << kShiftSector << ":" << ((mytr.sec & kMaskSector) << kShiftSector) << " SubSec " << mytr.subsec << ":" << kMaskSubSec << ":" << kShiftSubSec << ":" << ((mytr.subsec & kMaskSubSec) << kShiftSubSec) << " Index " << std::hex << indx << std::dec << std::endl; + + trformIndex_.push_back(indx); + trformTranX_.push_back(mytr.h3v.x()); + trformTranY_.push_back(mytr.h3v.y()); + trformTranZ_.push_back(mytr.h3v.z()); + trformRotXX_.push_back(mytr.hr.xx()); + trformRotYX_.push_back(mytr.hr.yx()); + trformRotZX_.push_back(mytr.hr.zx()); + trformRotXY_.push_back(mytr.hr.xy()); + trformRotYY_.push_back(mytr.hr.yy()); + trformRotZY_.push_back(mytr.hr.zy()); + trformRotXZ_.push_back(mytr.hr.xz()); + trformRotYZ_.push_back(mytr.hr.yz()); + trformRotZZ_.push_back(mytr.hr.zz()); +#ifdef DebugLog + unsigned int k = trformIndex_.size() - 1; + std::cout << "HGCalParameters[" << k << "] Index " << std::hex + << trformIndex_[k] << std::dec << " (" << mytr.zp << ", "<< mytr.lay + << ", " << mytr.sec << ", " << mytr.subsec << ") Translation (" + << trformTranX_[k] << ", " << trformTranY_[k] << ", " + << trformTranZ_[k] << ") Rotation (" << trformRotXX_[k] << ", " + << trformRotYX_[k] << ", " << trformRotZX_[k] << ", " + << trformRotXY_[k] << ", " << trformRotYY_[k] << ", " + << trformRotZY_[k] << ", " << trformRotXZ_[k] << ", " + << trformRotYZ_[k] << ", " << trformRotZZ_[k] << std::endl; +#endif +} + +HGCalParameters::hgtrform HGCalParameters::getTrForm(unsigned int k) const { + + HGCalParameters::hgtrform mytr; + if (k < trformIndex_.size()) { + int zp = ((trformIndex_[k] >> kShiftZside) & kMaskZside); + mytr.zp = (zp == 1) ? 1 : -1; + mytr.lay = ((trformIndex_[k] >> kShiftLayer) & kMaskLayer); + mytr.sec = ((trformIndex_[k] >> kShiftSector) & kMaskSector); + mytr.subsec= ((trformIndex_[k] >> kShiftSubSec) & kMaskSubSec); + mytr.h3v = CLHEP::Hep3Vector(trformTranX_[k],trformTranY_[k],trformTranZ_[k]); + const CLHEP::HepRep3x3 rotation(trformRotXX_[k],trformRotXY_[k],trformRotXZ_[k], + trformRotYX_[k],trformRotYY_[k],trformRotYZ_[k], + trformRotZX_[k],trformRotZY_[k],trformRotZZ_[k]); + mytr.hr = CLHEP::HepRotation(rotation); + } else { + mytr.zp = mytr.lay = mytr.sec = mytr.subsec = 0; + } +#ifdef DebugLog + std::cout << "HGCalParameters[" << k << "] Index " << std::hex + << trformIndex_[k] << std::dec << " (" << mytr.zp << ", "<< mytr.lay + << ", " << mytr.sec << ", " << mytr.subsec << ") Translation (" + << mytr.h3v.x() << ", " << mytr.h3v.y() << ", " << mytr.h3v.z() + << ") Rotation (" << mytr.hr.xx() << ", " << mytr.hr.yx() << ", " + << mytr.hr.zx() << ", " << mytr.hr.xy() << ", " << mytr.hr.yy() + << ", " << mytr.hr.zy() << ", " << mytr.hr.xz() << ", " + << mytr.hr.yz() << ", " << mytr.hr.zz() << std::endl; +#endif + return mytr; +} + +void HGCalParameters::addTrForm(const CLHEP::Hep3Vector& h3v) { + + unsigned int k = trformTranX_.size(); + if (k > 0) { + trformTranX_[k-1] += h3v.x(); + trformTranY_[k-1] += h3v.y(); + trformTranZ_[k-1] += h3v.z(); + } +} + +void HGCalParameters::scaleTrForm(double scale) { + unsigned int k = trformTranX_.size(); + if (k > 0) { + trformTranX_[k-1] *= scale; + trformTranY_[k-1] *= scale; + trformTranZ_[k-1] *= scale; + } +} diff --git a/CondFormats/GeometryObjects/src/classes_def.xml b/CondFormats/GeometryObjects/src/classes_def.xml index aaed2bca3a0cf..2bab4add55954 100644 --- a/CondFormats/GeometryObjects/src/classes_def.xml +++ b/CondFormats/GeometryObjects/src/classes_def.xml @@ -89,7 +89,7 @@ - - + diff --git a/Geometry/CaloGeometry/src/FlatTrd.cc b/Geometry/CaloGeometry/src/FlatTrd.cc index 75a7313ce15d8..7321d4a8bcddc 100644 --- a/Geometry/CaloGeometry/src/FlatTrd.cc +++ b/Geometry/CaloGeometry/src/FlatTrd.cc @@ -97,11 +97,21 @@ FlatTrd::~FlatTrd() {} GlobalPoint FlatTrd::getPosition(const Pt3D& local ) const { Pt3D glb = m_tr*local; +#ifdef DebugLog + std::cout << "FlatTrd::Local " << local.x() << ":" << local.y() << ":" + << local.z() << " Global " << glb.x() << ":" << glb.y() << ":" + << glb.z() << std::endl; +#endif return GlobalPoint(glb.x(),glb.y(),glb.z()); } Pt3D FlatTrd::getLocal(const GlobalPoint& global) const { Pt3D local = m_tr.inverse()*Pt3D(global.x(),global.y(),global.z()); +#ifdef DebugLog + std::cout << "FlatTrd::Global " << global.x() << ":" << global.y() << ":" + << global.z() << " Local " << local.x() << ":" << local.y() << ":" + << local.z() << std::endl; +#endif return local; } diff --git a/Geometry/HGCalCommonData/BuildFile.xml b/Geometry/HGCalCommonData/BuildFile.xml index 3fcdb8c74741a..eaa5971553706 100644 --- a/Geometry/HGCalCommonData/BuildFile.xml +++ b/Geometry/HGCalCommonData/BuildFile.xml @@ -1,5 +1,6 @@ + diff --git a/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h b/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h index 17f413022cd0a..78b383d8a09df 100644 --- a/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h +++ b/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h @@ -57,26 +57,25 @@ class HGCalDDDConstants { int numberCellsHexagon(int wafer) const; int sectors() const {return hgpar_->nSectors_;} std::pair simToReco(int cell, int layer, int mod, bool half) const; + unsigned int volumes() const {return hgpar_->moduleLayR_.size();} int waferFromCopy(int copy) const; bool waferInLayer(int wafer, int lay, bool reco) const; - GlobalPoint waferPosition(int wafer) const {return hgpar_->waferPos_.at(wafer); } + std::pair waferPosition(int wafer) const; int wafers() const; int waferToCopy(int wafer) const {return ((wafer>=0)&&(wafer< (int)(hgpar_->waferCopy_.size()))) ? hgpar_->waferCopy_[wafer] : (int)(hgpar_->waferCopy_.size());} int waferTypeT(int wafer) const {return ((wafer>=0)&&(wafer<(int)(hgpar_->waferTypeT_.size()))) ? hgpar_->waferTypeT_[wafer] : 0;} + HGCalParameters::hgtrap getModule(unsigned int k, bool hexType, bool reco) const; + std::vector getModules() const; - std::vector::const_iterator getFirstModule(bool reco=false) const {return (reco ? hgpar_->moduler_.begin() : hgpar_->modules_.begin());} - std::vector::const_iterator getLastModule(bool reco=false) const {return (reco ? hgpar_->moduler_.end() : hgpar_->modules_.end());} - std::vector::const_iterator getModule(int wafer) const; - std::vector::const_iterator getFirstTrForm() const {return hgpar_->trform_.begin();} - std::vector::const_iterator getLastTrForm() const {return hgpar_->trform_.end(); } - - const std::vector & getModules() const {return hgpar_->moduler_; } - const std::vector & getTrForms() const {return hgpar_->trform_; } + unsigned int getTrFormN() const {return hgpar_->trformIndex_.size();} + HGCalParameters::hgtrform getTrForm(unsigned int k) const {return hgpar_->getTrForm(k);} + std::vector getTrForms() const ; private: int cellHex(double xx, double yy, const double& cellR, - const std::vector& pos) const; + const std::vector& posX, + const std::vector& posY) const; std::pair getIndex(int lay, bool reco) const; void getParameterSquare(int lay, int subSec, bool reco, float& h, float& bl, float& tl, float& alpha) const; diff --git a/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.cc b/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.cc index e9f6a166e0897..20690337957b7 100644 --- a/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.cc +++ b/Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.cc @@ -13,7 +13,7 @@ #include "Geometry/HGCalCommonData/plugins/DDHGCalWaferAlgo.h" #include "CLHEP/Units/GlobalSystemOfUnits.h" -#define DebugLog +//#define DebugLog DDHGCalWaferAlgo::DDHGCalWaferAlgo() { #ifdef DebugLog diff --git a/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc b/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc index b75e8a3884167..dae8bd6cb8edb 100644 --- a/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc +++ b/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc @@ -116,15 +116,15 @@ std::pair HGCalDDDConstants::assignCellHexagon(float x, float y) const { double xx(x), yy(y); //First the wafer - int wafer = cellHex(xx, yy, rmax_, hgpar_->waferPos_); + int wafer = cellHex(xx, yy, rmax_, hgpar_->waferPosX_, hgpar_->waferPosY_); // Now the cell - xx -= hgpar_->waferPos_[wafer].x(); - yy -= hgpar_->waferPos_[wafer].y(); + xx -= hgpar_->waferPosX_[wafer]; + yy -= hgpar_->waferPosY_[wafer]; int cell(0); if (hgpar_->waferTypeT_[wafer] == 1) - cell = cellHex(xx, yy, 0.5*k_ScaleFromDDD*hgpar_->cellSize_[0], hgpar_->cellFine_); + cell = cellHex(xx, yy, 0.5*k_ScaleFromDDD*hgpar_->cellSize_[0], hgpar_->cellFineX_, hgpar_->cellFineY_); else - cell = cellHex(xx, yy, 0.5*k_ScaleFromDDD*hgpar_->cellSize_[1], hgpar_->cellCoarse_); + cell = cellHex(xx, yy, 0.5*k_ScaleFromDDD*hgpar_->cellSize_[1], hgpar_->cellCoarseX_, hgpar_->cellCoarseY_); return std::pair(wafer,cell); } @@ -181,18 +181,37 @@ std::pair HGCalDDDConstants::findCellSquare(int cell, float h, return std::pair(kx,ky); } -std::vector::const_iterator -HGCalDDDConstants::getModule(int wafer) const { - - std::vector::const_iterator itr; - int type = ((wafer>=0)&&(wafer<(int)(hgpar_->waferTypeL_.size()))) ? - hgpar_->waferTypeL_[wafer] : 3; - int indx(1); - for (itr = hgpar_->moduler_.begin(); itr != hgpar_->moduler_.end(); - ++itr, ++indx) { - if (type == indx) return itr; +HGCalParameters::hgtrap HGCalDDDConstants::getModule(unsigned int indx, + bool hexType, + bool reco) const { + + HGCalParameters::hgtrap mytr; + if (hexType) { + unsigned int type = (indx < hgpar_->waferTypeL_.size()) ? + hgpar_->waferTypeL_[indx] : 3; + if (type > 0) --type; + mytr = hgpar_->getModule(type, true); + } else { + if (reco) mytr = hgpar_->getModule(indx,true); + else mytr = hgpar_->getModule(indx,false); } - return hgpar_->moduler_.end(); + return mytr; +} + +std::vector HGCalDDDConstants::getModules() const { + + std::vector mytrs; + for (unsigned int k=0; kmoduleLayR_.size(); ++k) + mytrs.push_back(hgpar_->getModule(k,true)); + return mytrs; +} + +std::vector HGCalDDDConstants::getTrForms() const { + + std::vector mytrs; + for (unsigned int k=0; ktrformIndex_.size(); ++k) + mytrs.push_back(hgpar_->getTrForm(k)); + return mytrs; } bool HGCalDDDConstants::isValid(int lay, int mod, int cell, bool reco) const { @@ -211,7 +230,7 @@ bool HGCalDDDConstants::isValid(int lay, int mod, int cell, bool reco) const { (mod > 0 && mod <= modmax)); if (ok) { cellmax = (hgpar_->waferTypeT_[mod]==1) ? - (int)(hgpar_->cellFine_.size()) : (int)(hgpar_->cellCoarse_.size()); + (int)(hgpar_->cellFineX_.size()) : (int)(hgpar_->cellCoarseX_.size()); ok = (cell >=0 && cell <= cellmax); } } @@ -244,14 +263,14 @@ std::pair HGCalDDDConstants::locateCell(int cell, int lay, if (type != 1) x = -x; y = ((kxy.second+0.5)*cellSize-h); } else { - x = hgpar_->waferPos_[type].x(); - y = hgpar_->waferPos_[type].y(); + x = hgpar_->waferPosX_[type]; + y = hgpar_->waferPosY_[type]; if (hgpar_->waferTypeT_[type] == 1) { - x += hgpar_->cellFine_[cell].x(); - y += hgpar_->cellFine_[cell].y(); + x += hgpar_->cellFineX_[cell]; + y += hgpar_->cellFineY_[cell]; } else { - x += hgpar_->cellCoarse_[cell].x(); - y += hgpar_->cellCoarse_[cell].y(); + x += hgpar_->cellCoarseX_[cell]; + y += hgpar_->cellCoarseY_[cell]; } if (!reco) { x /= k_ScaleFromDDD; @@ -265,11 +284,11 @@ std::pair HGCalDDDConstants::locateCellHex(int cell, int wafer, bool reco) const { float x(0), y(0); if (hgpar_->waferTypeT_[wafer] == 1) { - x = hgpar_->cellFine_[cell].x(); - y = hgpar_->cellFine_[cell].y(); + x = hgpar_->cellFineX_[cell]; + y = hgpar_->cellFineY_[cell]; } else { - x = hgpar_->cellCoarse_[cell].x(); - y = hgpar_->cellCoarse_[cell].y(); + x = hgpar_->cellCoarseX_[cell]; + y = hgpar_->cellCoarseY_[cell]; } if (!reco) { x /= k_ScaleFromDDD; @@ -299,10 +318,10 @@ int HGCalDDDConstants::maxCells(int lay, bool reco) const { return maxCellsSquare(h, bl, tl, alpha, index.second); } else { unsigned int cells(0); - for (unsigned int k=0; kwaferPos_.size(); ++k) { + for (unsigned int k=0; kwaferTypeT_.size(); ++k) { if (waferInLayer(k,index.first)) { unsigned int cell = (hgpar_->waferTypeT_[k]==1) ? - (hgpar_->cellFine_.size()) : (hgpar_->cellCoarse_.size()); + (hgpar_->cellFineX_.size()) : (hgpar_->cellCoarseX_.size()); if (cell > cells) cells = cell; } } @@ -334,10 +353,10 @@ int HGCalDDDConstants::maxRows(int lay, bool reco) const { int i = index.first; if (i < 0) return kymax; if (geomMode() == HGCalGeometryMode::Square) { - float h = (reco) ? hgpar_->moduler_[i].h : hgpar_->modules_[i].h; + float h = (reco) ? hgpar_->moduleHR_[i] : hgpar_->moduleHS_[i]; kymax = floor((2*h)/index.second); } else { - for (unsigned int k=0; kwaferPos_.size(); ++k) { + for (unsigned int k=0; kwaferCopy_.size(); ++k) { if (waferInLayer(k,i)) { int ky = ((hgpar_->waferCopy_[k])/100)%100; if (ky > kymax) kymax = ky; @@ -351,7 +370,7 @@ int HGCalDDDConstants::modules(int lay, bool reco) const { int nmod(0); std::pair index = getIndex(lay, reco); if (index.first < 0) return nmod; - for (unsigned int k=0; kwaferPos_.size(); ++k) { + for (unsigned int k=0; kwaferPosX_.size(); ++k) { if (waferInLayer(k,index.first)) ++nmod; } return nmod; @@ -399,13 +418,13 @@ int HGCalDDDConstants::newCell(int kx, int ky, int lay, int subSec) const { std::pair index = getIndex(lay, true); int i = index.first; if (i < 0) return maxCells(true); - float alpha = (subSec == 0) ? hgpar_->modules_[i].alpha : subSec; + float alpha = (subSec == 0) ? hgpar_->moduleAlphaS_[i] : subSec; float cellSize = index.second; float a = (alpha==0) ? - (2*hgpar_->moduler_[i].h/(hgpar_->moduler_[i].tl-hgpar_->moduler_[i].bl)) : - (hgpar_->moduler_[i].h/(hgpar_->moduler_[i].tl-hgpar_->moduler_[i].bl)); - float b = 2*hgpar_->moduler_[i].h*hgpar_->moduler_[i].bl/ - (hgpar_->moduler_[i].tl-hgpar_->moduler_[i].bl); + (2*hgpar_->moduleHR_[i]/(hgpar_->moduleTlR_[i]-hgpar_->moduleBlR_[i])) : + (hgpar_->moduleHR_[i]/(hgpar_->moduleTlR_[i]-hgpar_->moduleBlR_[i])); + float b = 2*hgpar_->moduleHR_[i]*hgpar_->moduleBlR_[i]/ + (hgpar_->moduleTlR_[i]-hgpar_->moduleBlR_[i]); int icell(kx); for (int iky=0; iky HGCalDDDConstants::numberCells(int lay, bool reco) const { getParameterSquare(i,0,reco,h,bl,tl,alpha); return numberCellsSquare(h, bl, tl, alpha, index.second); } else { - for (unsigned int k=0; kwaferPos_.size(); ++k) { + for (unsigned int k=0; kwaferTypeT_.size(); ++k) { if (waferInLayer(k,i)) { unsigned int cell = (hgpar_->waferTypeT_[k]==1) ? - (hgpar_->cellFine_.size()) : (hgpar_->cellCoarse_.size()); + (hgpar_->cellFineX_.size()) : (hgpar_->cellCoarseX_.size()); ncell.push_back((int)(cell)); } } @@ -453,9 +472,9 @@ int HGCalDDDConstants::numberCellsHexagon(int wafer) const { int ncell(0); if (wafer >= 0 && wafer < (int)(hgpar_->waferTypeT_.size())) { if (hgpar_->waferTypeT_[wafer]==1) - ncell = (int)(hgpar_->cellFine_.size()); + ncell = (int)(hgpar_->cellFineX_.size()); else - ncell = (int)(hgpar_->cellCoarse_.size()); + ncell = (int)(hgpar_->cellCoarseX_.size()); } return ncell; } @@ -468,12 +487,12 @@ std::pair HGCalDDDConstants::simToReco(int cell, int lay, int mod, if (i < 0) return std::pair(-1,-1); int kx(-1), depth(-1); if (geomMode() == HGCalGeometryMode::Square) { - float h = hgpar_->modules_[i].h; - float bl = hgpar_->modules_[i].bl; - float tl = hgpar_->modules_[i].tl; + float h = hgpar_->moduleHS_[i]; + float bl = hgpar_->moduleBlS_[i]; + float tl = hgpar_->moduleTlS_[i]; float cellSize = hgpar_->cellFactor_[i]*index.second; - std::pair kxy = findCellSquare(cell, h, bl, tl, hgpar_->modules_[i].alpha, index.second); + std::pair kxy = findCellSquare(cell, h, bl, tl, hgpar_->moduleAlphaS_[i], index.second); depth = hgpar_->layerGroup_[i]; if (depth<0) return std::pair(-1,-1); kx = kxy.first/hgpar_->cellFactor_[i]; @@ -521,6 +540,18 @@ bool HGCalDDDConstants::waferInLayer(int wafer, int lay, bool reco) const { return waferInLayer(wafer,indx.first); } +std::pair HGCalDDDConstants::waferPosition(int wafer) const { + + std::pair xy; + if (wafer >= 0 && wafer < (int)(hgpar_->waferPosX_.size())) { + xy = std::pair(hgpar_->waferPosX_[wafer],hgpar_->waferPosY_[wafer]); + } else { + xy = std::pair(0,0); + } + return xy; +} + + int HGCalDDDConstants::wafers() const { int wafer(0); @@ -532,13 +563,14 @@ int HGCalDDDConstants::wafers() const { } int HGCalDDDConstants::cellHex(double xx, double yy, const double& cellR, - const std::vector& pos) const { + const std::vector& posX, + const std::vector& posY) const { int num(0); const double tol(0.00001); double cellY = 2.0*cellR*tan30deg_; - for (unsigned int k=0; k HGCalDDDConstants::getIndex(int lay, bool reco) const { float cell(0); if (geomMode() == HGCalGeometryMode::Square) { indx = (reco ? hgpar_->depthIndex_[lay-1] : hgpar_->layerIndex_[lay-1]); - cell = (reco ? hgpar_->moduler_[indx].cellSize : hgpar_->modules_[indx].cellSize); + cell = (reco ? hgpar_->moduleCellR_[indx] : hgpar_->moduleCellS_[indx]); } else { indx = (reco ? hgpar_->depthLayerF_[lay-1] : hgpar_->layerIndex_[lay-1]); - cell = (reco ? hgpar_->moduler_[0].cellSize : hgpar_->modules_[0].cellSize); + cell = (reco ? hgpar_->moduleCellR_[0] : hgpar_->moduleCellS_[0]); } return std::pair(indx,cell); } @@ -570,23 +602,24 @@ void HGCalDDDConstants::getParameterSquare(int lay, int subSec, bool reco, float& h, float& bl, float& tl, float& alpha) const { if (reco) { - h = hgpar_->moduler_[lay].h; - bl = hgpar_->moduler_[lay].bl; - tl = hgpar_->moduler_[lay].tl; - alpha= hgpar_->moduler_[lay].alpha; + h = hgpar_->moduleHR_[lay]; + bl = hgpar_->moduleBlR_[lay]; + tl = hgpar_->moduleTlR_[lay]; + alpha= hgpar_->moduleAlphaR_[lay]; if ((subSec>0 && alpha<0) || (subSec<=0 && alpha>0)) alpha = -alpha; } else { - h = hgpar_->modules_[lay].h; - bl = hgpar_->modules_[lay].bl; - tl = hgpar_->modules_[lay].tl; - alpha= hgpar_->modules_[lay].alpha; + h = hgpar_->moduleHS_[lay]; + bl = hgpar_->moduleBlS_[lay]; + tl = hgpar_->moduleTlS_[lay]; + alpha= hgpar_->moduleAlphaS_[lay]; } } bool HGCalDDDConstants::waferInLayer(int wafer, int lay) const { double rr = 2*rmax_*tan30deg_; - double rpos = hgpar_->waferPos_[wafer].perp(); + double rpos = std::sqrt(hgpar_->waferPosX_[wafer]*hgpar_->waferPosX_[wafer]+ + hgpar_->waferPosY_[wafer]*hgpar_->waferPosY_[wafer]); bool in = (rpos-rr >= hgpar_->rMinLayHex_[lay] && rpos+rr <= hgpar_->rMaxLayHex_[lay]); return in; diff --git a/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc b/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc index 56e8c51cccea3..f90624e6e9a95 100644 --- a/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc +++ b/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc @@ -39,6 +39,7 @@ void HGCalGeomParameters::loadGeometrySquare(const DDFilteredView& _fv, bool dodet(true), first(true); int zpFirst(0); std::vector trforms; + std::vector trformUse; while (dodet) { const DDSolid & sol = fv.logicalPart().solid(); @@ -66,7 +67,7 @@ void HGCalGeomParameters::loadGeometrySquare(const DDFilteredView& _fv, break; } } - php.modules_.push_back(mytr); + php.fillModule(mytr,false); if (php.layer_.size() == 0) php.nSectors_ = 1; php.layer_.push_back(lay); } else if (std::find(php.layer_.begin(),php.layer_.end(),lay) == @@ -89,8 +90,8 @@ void HGCalGeomParameters::loadGeometrySquare(const DDFilteredView& _fv, mytrf.subsec= subs; mytrf.h3v = h3v; mytrf.hr = hr; - mytrf.used = false; trforms.push_back(mytrf); + trformUse.push_back(false); } dodet = fv.next(); } @@ -116,11 +117,10 @@ void HGCalGeomParameters::loadGeometrySquare(const DDFilteredView& _fv, for (unsigned int i=0; i 0) { - php.trform_.back().h3v /= nz; + php.scaleTrForm(double(1.0/nz)); } } } } #ifdef DebugLog - std::cout << "Obtained " << php.trform_.size() << " transformation matrices" - << std::endl; - for (unsigned int k=0; k layers; std::vector trforms; + std::vector trformUse; while (dodet) { const DDSolid & sol = fv.logicalPart().solid(); @@ -258,8 +264,8 @@ void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv, mytrf.subsec= 0; mytrf.h3v = h3v; mytrf.hr = hr; - mytrf.used = false; trforms.push_back(mytrf); + trformUse.push_back(false); } } dodet = fv.next(); @@ -316,7 +322,7 @@ void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv, mytr.tl = php.waferR_; mytr.h = php.waferR_; mytr.dz = dz; mytr.alpha = 0.0; mytr.cellSize = waferSize_; - php.modules_.push_back(mytr); + php.fillModule(mytr,false); names.push_back(name); } } @@ -413,23 +419,23 @@ void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv, } for (unsigned int i=0; i 0) { - php.trform_.back().h3v /= nz; + php.scaleTrForm(double(1.0/nz)); } } } @@ -439,7 +445,8 @@ void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv, for (std::map::iterator itr = wafers.begin(); itr != wafers.end(); ++itr) { php.waferCopy_.push_back(itr->first); - php.waferPos_.push_back(itr->second); + php.waferPosX_.push_back(itr->second.x()); + php.waferPosY_.push_back(itr->second.y()); std::map::iterator ktr = wafertype.find(itr->first); int typet = (ktr == wafertype.end()) ? 0 : (ktr->second); php.waferTypeT_.push_back(typet); @@ -470,7 +477,8 @@ void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv, if (xy.first > 0) xy.first -= 0.001; else xy.first += 0.001; } - php.cellFine_.push_back(GlobalPoint(xy.first,xy.second,0)); + php.cellFineX_.push_back(xy.first); + php.cellFineY_.push_back(xy.second); } } itrf = wafers.end(); @@ -489,7 +497,8 @@ void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv, if (xy.first > 0) xy.first -= 0.001; else xy.first += 0.001; } - php.cellCoarse_.push_back(GlobalPoint(xy.first,xy.second,0)); + php.cellCoarseX_.push_back(xy.first); + php.cellCoarseY_.push_back(xy.second); } } int depth(0); @@ -507,16 +516,17 @@ void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv, } } } - php.moduler_.push_back(php.modules_[0]); - php.moduler_.back().bl *= k_ScaleFromDDD; - php.moduler_.back().tl *= k_ScaleFromDDD; - php.moduler_.back().h *= k_ScaleFromDDD; - php.moduler_.back().dz *= k_ScaleFromDDD; - double dz = php.moduler_.back().dz; - php.moduler_.push_back(php.moduler_[0]); - php.moduler_.back().dz = 2*dz; - php.moduler_.push_back(php.moduler_[0]); - php.moduler_.back().dz = 3*dz; + HGCalParameters::hgtrap mytr = php.getModule(0, false); + mytr.bl *= k_ScaleFromDDD; + mytr.tl *= k_ScaleFromDDD; + mytr.h *= k_ScaleFromDDD; + mytr.dz *= k_ScaleFromDDD; + double dz = mytr.dz; + php.fillModule(mytr, true); + mytr.dz = 2*dz; + php.fillModule(mytr, true); + mytr.dz = 3*dz; + php.fillModule(mytr, true); #ifdef DebugLog std::cout << "HGCalGeomParameters: rmax = " << rmax << ", ymax = " << ymax << std::endl; @@ -539,38 +549,41 @@ void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv, for (unsigned int i=0; i xy; @@ -145,15 +145,14 @@ GlobalPoint HGCalGeometry::getPosition(const DetId& id) const { xy = topology().dddConstants().locateCellHex(id_.iCell,id_.iSec,true); } const HepGeom::Point3D lcoord(xy.first,xy.second,0); + glob = m_cellVec[cellIndex].getPosition(lcoord); #ifdef DebugLog - std::cout << "getPosition:: index " << cellIndex << " Local " << xy.first - << ":" << xy.second << " ID " << id_.iCell << ":" << id_.iLay - << " Global " << m_cellVec[cellIndex].getPosition(lcoord) - << " Cell" << m_cellVec[cellIndex]; + std::cout << "getPosition:: index " << cellIndex << " Local " << lcoord.x() + << ":" << lcoord.y() << " ID " << id_.iCell << ":" << id_.iLay + << " Global " << glob << std::endl; #endif - return m_cellVec[cellIndex].getPosition(lcoord); } - return GlobalPoint(); + return glob; } HGCalGeometry::CornersVec HGCalGeometry::getCorners(const DetId& id) const { @@ -206,7 +205,7 @@ DetId HGCalGeometry::getClosestCell(const GlobalPoint& r) const { #ifdef DebugLog std::cout << "getClosestCell: local " << local << " Id " << id_.zside << ":" << id_.iLay << ":" << id_.iSec << ":" << id_.iSubSec - << ":" << id_.iCell << " Cell " << m_cellVec[cellIndex]; + << ":" << id_.iCell << std::endl; #endif //check if returned cell is valid @@ -257,7 +256,7 @@ const CaloCellGeometry* HGCalGeometry::cellGeomPtr(uint32_t index) const { return 0; const CaloCellGeometry* cell ( &m_cellVec[ index ] ) ; #ifdef DebugLog - std::cout << "cellGeomPtr " << m_cellVec[index]; + // std::cout << "cellGeomPtr " << m_cellVec[index]; #endif if (0 == cell->param()) return 0; return cell; @@ -316,9 +315,7 @@ namespace }; } -void -HGCalGeometry::sortDetIds( void ) -{ +void HGCalGeometry::sortDetIds( void ) { m_validIds.shrink_to_fit(); std::sort( m_validIds.begin(), m_validIds.end(), rawIdSort()); } @@ -327,9 +324,9 @@ void HGCalGeometry::getSummary( CaloSubdetectorGeometry::TrVec& trVector, CaloSubdetectorGeometry::IVec& iVector, CaloSubdetectorGeometry::DimVec& dimVector, - CaloSubdetectorGeometry::IVec& dinsVector ) const -{ - unsigned int numberOfCells = m_topology.totalGeomModules(); // total Geom Modules both sides + CaloSubdetectorGeometry::IVec& dinsVector ) const { + + unsigned int numberOfCells = topology().totalGeomModules(); // total Geom Modules both sides unsigned int numberOfShapes = HGCalGeometry::k_NumberOfShapes; unsigned int numberOfParametersPerShape = HGCalGeometry::k_NumberOfParametersPerShape; @@ -338,19 +335,38 @@ HGCalGeometry::getSummary( CaloSubdetectorGeometry::TrVec& trVector, dimVector.reserve( numberOfShapes * numberOfParametersPerShape ); dinsVector.reserve( numberOfCells ); - for( auto volItr = m_topology.dddConstants().getFirstModule( true ); - volItr != m_topology.dddConstants().getLastModule( true ); ++volItr ) - { - ParmVec params( HGCalGeometry::k_NumberOfParametersPerShape, 0 ); - params[0] = volItr->dz; - params[1] = params[2] = 0; - params[3] = params[7] = volItr->h; - params[4] = params[8] = volItr->bl; - params[5] = params[9] = volItr->tl; - params[6] = params[10]= volItr->alpha; - params[11]= volItr->cellSize; - - dimVector.insert( dimVector.end(), params.begin(), params.end()); + if (topology().geomMode() == HGCalGeometryMode::Square) { + for (unsigned int k=0; k ::const_iterator trItr; - std::vector::const_iterator volItr; ParmVec params(HGCalGeometry::k_NumberOfParametersPerShape,0); unsigned int counter(0); - for (trItr = topology.dddConstants().getFirstTrForm(); - trItr != topology.dddConstants().getLastTrForm(); ++trItr) { - int zside = trItr->zp; - int layer = trItr->lay; +#ifdef DebugLog + std::cout << "HGCalGeometryLoader with # of transformation matrices " + << topology.dddConstants().getTrFormN() << " and " + << topology.dddConstants().volumes() << ":" + << topology.dddConstants().sectors() << " volumes" << std::endl; +#endif + for (unsigned itr=0; itrsec; - int subSec = (detType ? trItr->subsec : 0); - const HepGeom::Transform3D ht3d (trItr->hr, trItr->h3v); + int sector = mytr.sec; + int subSec = (detType ? mytr.subsec : 0); + const HepGeom::Transform3D ht3d (mytr.hr, mytr.h3v); DetId detId= ((subdet == HGCEE) ? (DetId)(HGCEEDetId(subdet,zside,layer,sector,subSec,0)) : (DetId)(HGCHEDetId(subdet,zside,layer,sector,subSec,0))); @@ -58,19 +62,18 @@ HGCalGeometry* HGCalGeometryLoader::build (const HGCalTopology& topology) { << subSec << " transf " << ht3d.getTranslation() << " and " << ht3d.getRotation(); #endif - for (volItr = topology.dddConstants().getFirstModule(true); - volItr != topology.dddConstants().getLastModule(true); ++volItr) { - if (volItr->lay == layer) { - double alpha = ((detType && subSec == 0) ? -fabs(volItr->alpha) : - fabs(volItr->alpha)); - params[0] = volItr->dz; + for (unsigned int k=0; kh; - params[4] = params[8] = volItr->bl; - params[5] = params[9] = volItr->tl; + params[3] = params[7] = vol.h; + params[4] = params[8] = vol.bl; + params[5] = params[9] = vol.tl; params[6] = params[10]= alpha; - params[11]= volItr->cellSize; - + params[11]= vol.cellSize; buildGeom(params, ht3d, detId, geom); counter++; break; @@ -82,21 +85,21 @@ HGCalGeometry* HGCalGeometryLoader::build (const HGCalTopology& topology) { int type = topology.dddConstants().waferTypeT(wafer); if (type != 1) type = 0; DetId detId = (DetId)(HGCalDetId(subdet,zside,layer,type,wafer,0)); - GlobalPoint w = topology.dddConstants().waferPosition(wafer); - double xx = (zside > 0) ? w.x() : -w.x(); - CLHEP::Hep3Vector h3v(xx,w.y(),trItr->h3v.z()); - const HepGeom::Transform3D ht3d (trItr->hr, h3v); + std::pair w = topology.dddConstants().waferPosition(wafer); + double xx = (zside > 0) ? w.first : -w.first; + CLHEP::Hep3Vector h3v(xx,w.second,mytr.h3v.z()); + const HepGeom::Transform3D ht3d (mytr.hr, h3v); #ifdef DebugLog std::cout << "HGCalGeometryLoader:: Wafer:Type " << wafer << ":" << type << " transf " << ht3d.getTranslation() << " and " << ht3d.getRotation(); #endif - volItr = topology.dddConstants().getModule(wafer); - params[0] = volItr->dz; + HGCalParameters::hgtrap vol = topology.dddConstants().getModule(wafer,true,true); + params[0] = vol.dz; params[1] = params[2] = 0; - params[3] = params[7] = volItr->h; - params[4] = params[8] = volItr->bl; - params[5] = params[9] = volItr->tl; + params[3] = params[7] = vol.h; + params[4] = params[8] = vol.bl; + params[5] = params[9] = vol.tl; params[6] = params[10]= 0; params[11]= topology.dddConstants().cellSizeHex(type); diff --git a/Geometry/HGCalGeometry/test/HGCalGeometryTester.cc b/Geometry/HGCalGeometry/test/HGCalGeometryTester.cc index 56f85ba32e4c9..088ad7e796df5 100644 --- a/Geometry/HGCalGeometry/test/HGCalGeometryTester.cc +++ b/Geometry/HGCalGeometry/test/HGCalGeometryTester.cc @@ -71,13 +71,13 @@ void HGCalGeometryTester::doTest(const HGCalGeometry& geom, int zsides[] = {1, -1}; int cells[] = {1, 51, 101}; int wafers[] = {1, 101, 201, 301, 401}; - int types[] = {1, 0, 1, 0, 1}; int ismax = (squareCell) ? 3 : 5; for (int iz = 0; iz < 2; ++iz) { int zside = zsides[iz]; for (int is = 0; is < ismax; ++is) { int sector = (squareCell) ? sectors[is] : wafers[is]; - int type = (squareCell) ? 0 : types[is]; + int type = (squareCell) ? 0 : geom.topology().dddConstants().waferTypeT(sector); + if (type != 1) type = 0; for (int il = 0; il < 3; ++il) { int layer = layers[il]; for (int ic = 0; ic < 3; ++ic) { diff --git a/SimDataFormats/CaloTest/src/HGCalTestNumbering.cc b/SimDataFormats/CaloTest/src/HGCalTestNumbering.cc index 69c245f809424..79b95d60a60bf 100644 --- a/SimDataFormats/CaloTest/src/HGCalTestNumbering.cc +++ b/SimDataFormats/CaloTest/src/HGCalTestNumbering.cc @@ -1,7 +1,7 @@ #include "SimDataFormats/CaloTest/interface/HGCalTestNumbering.h" #include -#define DebugLog +//#define DebugLog uint32_t HGCalTestNumbering::packSquareIndex(int zp, int lay, int sec, int subsec, int cell) { From aafd23bab0897738752ff44709381b0f57108417 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Sun, 27 Dec 2015 19:18:14 +0100 Subject: [PATCH 3/6] Modify the cff's with HGCal --- .../python/GeometryExtended2023Reco_cff.py | 25 ++++++++++++++++--- .../python/GeometryExtended2023_cff.py | 1 + 2 files changed, 23 insertions(+), 3 deletions(-) diff --git a/Configuration/Geometry/python/GeometryExtended2023Reco_cff.py b/Configuration/Geometry/python/GeometryExtended2023Reco_cff.py index 1f8c24b48bc4a..10cf5f8f3ac94 100644 --- a/Configuration/Geometry/python/GeometryExtended2023Reco_cff.py +++ b/Configuration/Geometry/python/GeometryExtended2023Reco_cff.py @@ -17,6 +17,7 @@ from Geometry.MuonNumbering.muonNumberingInitialization_cfi import * from RecoMuon.DetLayers.muonDetLayerGeometry_cfi import * from Geometry.GEMGeometryBuilder.gemGeometry_cfi import * +from Geometry.GEMGeometryBuilder.me0Geometry_cfi import * # Alignment from Geometry.TrackerGeometryBuilder.idealForDigiTrackerGeometry_cff import * @@ -25,10 +26,28 @@ trackerGeometry.applyAlignment = cms.bool(False) # Calorimeters +from Geometry.CaloEventSetup.HGCalTopology_cfi import * +from Geometry.HGCalGeometry.HGCalGeometryESProducer_cfi import * from Geometry.CaloEventSetup.CaloTopology_cfi import * -from Geometry.CaloEventSetup.CaloGeometry_cff import * +from Geometry.CaloEventSetup.CaloGeometryBuilder_cfi import * + +CaloGeometryBuilder = cms.ESProducer("CaloGeometryBuilder", + SelectedCalos = cms.vstring('HCAL' , + 'ZDC' , + 'CASTOR' , + 'EcalBarrel' , + 'TOWER' ) +) + +from Geometry.EcalAlgo.EcalBarrelGeometry_cfi import * +from Geometry.HcalEventSetup.HcalGeometry_cfi import * +from Geometry.HcalEventSetup.CaloTowerGeometry_cfi import * +from Geometry.HcalEventSetup.CaloTowerTopology_cfi import * +from Geometry.HcalCommonData.hcalDDDRecConstants_cfi import * +from Geometry.HcalEventSetup.hcalTopologyIdeal_cfi import * +from Geometry.ForwardGeometry.ForwardGeometry_cfi import * + from Geometry.CaloEventSetup.EcalTrigTowerConstituents_cfi import * from Geometry.EcalMapping.EcalMapping_cfi import * from Geometry.EcalMapping.EcalMappingRecord_cfi import * -from Geometry.HcalCommonData.hcalDDDRecConstants_cfi import * -from Geometry.HcalEventSetup.hcalTopologyIdeal_cfi import * + diff --git a/Configuration/Geometry/python/GeometryExtended2023_cff.py b/Configuration/Geometry/python/GeometryExtended2023_cff.py index 3bf61ea7fe648..d494419f11c1c 100644 --- a/Configuration/Geometry/python/GeometryExtended2023_cff.py +++ b/Configuration/Geometry/python/GeometryExtended2023_cff.py @@ -8,5 +8,6 @@ from Geometry.TrackerNumberingBuilder.trackerNumberingGeometry_cfi import * from Geometry.HcalCommonData.hcalParameters_cfi import * from Geometry.HcalCommonData.hcalDDDSimConstants_cfi import * +from Geometry.HGCalCommonData.hgcalParametersInitialization_cfi import * from Geometry.HGCalCommonData.hgcalNumberingInitialization_cfi import * from Geometry.CaloEventSetup.HGCalTopology_cfi import * From 9ef390e95eac6fad97e55b2b6f4e9f03fae0e944 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Mon, 28 Dec 2015 15:44:26 +0100 Subject: [PATCH 4/6] Update the geometry file usage in test --- .../HGCalRecProducers/test/testHGCalRecoLocal_cfg.py | 4 ++-- SimCalorimetry/HGCalSimProducers/test/testHGCalDigi_cfg.py | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/RecoLocalCalo/HGCalRecProducers/test/testHGCalRecoLocal_cfg.py b/RecoLocalCalo/HGCalRecProducers/test/testHGCalRecoLocal_cfg.py index ce21827f1b70a..5d6acbebbcd0d 100644 --- a/RecoLocalCalo/HGCalRecProducers/test/testHGCalRecoLocal_cfg.py +++ b/RecoLocalCalo/HGCalRecProducers/test/testHGCalRecoLocal_cfg.py @@ -11,8 +11,8 @@ process.load('FWCore.MessageService.MessageLogger_cfi') process.load('Configuration.EventContent.EventContent_cff') process.load('SimGeneral.MixingModule.mixNoPU_cfi') -process.load('Configuration.Geometry.GeometryExtended2023DevReco_cff') -process.load('Configuration.Geometry.GeometryExtended2023Dev_cff') +process.load('Configuration.Geometry.GeometryExtended2023Reco_cff') +process.load('Configuration.Geometry.GeometryExtended2023_cff') process.load('Configuration.StandardSequences.MagneticField_38T_PostLS1_cff') process.load('Configuration.StandardSequences.Generator_cff') process.load('IOMC.EventVertexGenerators.VtxSmearedGauss_cfi') diff --git a/SimCalorimetry/HGCalSimProducers/test/testHGCalDigi_cfg.py b/SimCalorimetry/HGCalSimProducers/test/testHGCalDigi_cfg.py index 0f2761d0d7ed0..22f8f9cdc3648 100644 --- a/SimCalorimetry/HGCalSimProducers/test/testHGCalDigi_cfg.py +++ b/SimCalorimetry/HGCalSimProducers/test/testHGCalDigi_cfg.py @@ -11,8 +11,8 @@ process.load('FWCore.MessageService.MessageLogger_cfi') process.load('Configuration.EventContent.EventContent_cff') process.load('SimGeneral.MixingModule.mixNoPU_cfi') -process.load('Configuration.Geometry.GeometryExtended2023DevReco_cff') -process.load('Configuration.Geometry.GeometryExtended2023Dev_cff') +process.load('Configuration.Geometry.GeometryExtended2023Reco_cff') +process.load('Configuration.Geometry.GeometryExtended2023_cff') process.load('Configuration.StandardSequences.MagneticField_38T_PostLS1_cff') process.load('Configuration.StandardSequences.Generator_cff') process.load('IOMC.EventVertexGenerators.VtxSmearedGauss_cfi') From 2be0fda003e42e5de5c2d21dd9376713aa87b069 Mon Sep 17 00:00:00 2001 From: Sunanda Date: Wed, 30 Dec 2015 20:22:23 +0100 Subject: [PATCH 5/6] Bug fix --- Geometry/HGCalCommonData/src/HGCalGeomParameters.cc | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) diff --git a/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc b/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc index f90624e6e9a95..1917f3dafa372 100644 --- a/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc +++ b/Geometry/HGCalCommonData/src/HGCalGeomParameters.cc @@ -409,7 +409,7 @@ void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv, for (std::map::iterator itr = layers.begin(); itr != layers.end(); ++itr) { if (itr->first == (int)(i+1)) { - php.layerIndex_.push_back(itr->first); + php.layerIndex_.push_back(i); php.rMinLayHex_.push_back(itr->second.rmin); php.rMaxLayHex_.push_back(itr->second.rmax); php.zLayerHex_.push_back(itr->second.zpos); @@ -532,10 +532,13 @@ void HGCalGeomParameters::loadGeometryHexagon(const DDFilteredView& _fv, << std::endl; std::cout << "HGCalGeomParameters finds " << php.zLayerHex_.size() << " layers" << std::endl; - for (unsigned int i=0; i Date: Wed, 30 Dec 2015 21:19:23 +0100 Subject: [PATCH 6/6] Needed for validation work --- Geometry/HGCalCommonData/interface/HGCalDDDConstants.h | 1 + Geometry/HGCalCommonData/src/HGCalDDDConstants.cc | 6 ++++++ 2 files changed, 7 insertions(+) diff --git a/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h b/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h index 78b383d8a09df..29d85cafb6188 100644 --- a/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h +++ b/Geometry/HGCalCommonData/interface/HGCalDDDConstants.h @@ -64,6 +64,7 @@ class HGCalDDDConstants { int wafers() const; int waferToCopy(int wafer) const {return ((wafer>=0)&&(wafer< (int)(hgpar_->waferCopy_.size()))) ? hgpar_->waferCopy_[wafer] : (int)(hgpar_->waferCopy_.size());} int waferTypeT(int wafer) const {return ((wafer>=0)&&(wafer<(int)(hgpar_->waferTypeT_.size()))) ? hgpar_->waferTypeT_[wafer] : 0;} + double waferZ(int layer, bool reco) const; HGCalParameters::hgtrap getModule(unsigned int k, bool hexType, bool reco) const; std::vector getModules() const; diff --git a/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc b/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc index dae8bd6cb8edb..b4e18746ba2aa 100644 --- a/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc +++ b/Geometry/HGCalCommonData/src/HGCalDDDConstants.cc @@ -551,6 +551,12 @@ std::pair HGCalDDDConstants::waferPosition(int wafer) const { return xy; } +double HGCalDDDConstants::waferZ(int lay, bool reco) const { + std::pair index = getIndex(lay, reco); + int i = index.first; + if (i < 0) return 0; + else return hgpar_->zLayerHex_[i]; +} int HGCalDDDConstants::wafers() const {