17 #include "CLHEP/Units/GlobalPhysicalConstants.h"
18 #include "CLHEP/Units/GlobalSystemOfUnits.h"
19 #include <unordered_set>
27 std::cout <<
"HGCalGeomParameters::HGCalGeomParameters() constructor\n";
33 std::cout <<
"HGCalGeomParameters::destructed!!!\n";
42 bool dodet(
true),
first(
true);
44 std::vector<HGCalParameters::hgtrform> trforms;
45 std::vector<bool> trformUse;
50 int isd = (name.find(sdTag) == std::string::npos) ? -1 : 1;
53 int nsiz = (int)(copy.size());
54 int lay = (nsiz > 0) ? copy[nsiz-1] : -1;
55 int sec = (nsiz > 1) ? copy[nsiz-2] : -1;
56 int zp = (nsiz > 3) ? copy[nsiz-4] : -1;
58 if (first) {first =
false; zpFirst = zp;}
61 mytr.
lay = lay; mytr.
bl = trp.
x1(); mytr.
tl = trp.
x2();
64 int subs = (trp.
alpha1()>0 ? 1 : 0);
68 if (lay == (
int)(
k+1)) {
81 fv.
rotation().GetComponents( x, y, z ) ;
82 const CLHEP::HepRep3x3
rotation ( x.X(), y.X(), z.X(),
84 x.Z(), y.Z(), z.Z() );
96 trforms.push_back(mytrf);
97 trformUse.push_back(
false);
102 edm::LogError(
"HGCalGeom") <<
"HGCalGeomParameters : mismatch in # of bins "
104 <<
" between geometry and specpar";
105 throw cms::Exception(
"DDException") <<
"HGCalGeomParameters: mismatch between geometry and specpar";
107 for (
unsigned int i=0;
i<php.
layer_.size(); ++
i) {
108 for (
unsigned int k=0;
k<php.
layer_.size(); ++
k) {
117 <<
" modules for " << sdTag <<
" with " << php.
nSectors_
118 <<
" sectors and " << trforms.size() <<
" transformation matrices"
131 for (
unsigned int i=0;
i<php.
layer_.size(); ++
i) {
160 <<
" depths" << std::endl;
163 std::cout <<
"Module[" <<
i <<
":" << k <<
"] First Layer "
171 for (
unsigned int i=0;
i<php.
layer_.size(); ++
i) {
172 for (
unsigned int i1=0; i1<trforms.size(); ++i1) {
173 if (!trformUse[i1] && php.
layerGroup_[trforms[i1].lay-1] ==
176 trforms[i1].lay = (
i+1);
177 trformUse[i1] =
true;
180 for (
unsigned int i2=i1+1; i2<trforms.size(); ++i2) {
181 if (!trformUse[i2] && trforms[i2].zp == trforms[i1].zp &&
183 trforms[i2].sec == trforms[i1].sec &&
184 trforms[i2].subsec == trforms[i1].subsec) {
187 trformUse[i2] =
true;
198 <<
" transformation matrices" << std::endl;
221 std::map<int,HGCalGeomParameters::layerParameters>
layers;
222 std::vector<HGCalParameters::hgtrform> trforms;
223 std::vector<bool> trformUse;
230 int nsiz = (int)(copy.size());
231 int lay = (nsiz > 0) ? copy[nsiz-1] : 0;
232 int zp = (nsiz > 2) ? copy[nsiz-3] : -1;
233 if (zp != 1) zp = -1;
235 edm::LogError(
"HGCalGeom") <<
"Funny layer # " << lay <<
" zp "
236 << zp <<
" in " << nsiz <<
" components";
241 std::map<int,HGCalGeomParameters::layerParameters>::iterator itr = layers.find(lay);
242 if (itr == layers.end()) {
248 layers[lay] = laypar;
251 fv.
rotation().GetComponents( x, y, z ) ;
252 const CLHEP::HepRep3x3
rotation ( x.X(), y.X(), z.X(),
254 x.Z(), y.Z(), z.Z() );
260 const CLHEP::Hep3Vector h3v ( xx, yy, fv.
translation().Z() );
268 trforms.push_back(mytrf);
269 trformUse.push_back(
false);
278 std::unordered_map<int32_t,int32_t> copies;
280 std::vector<int32_t> wafer2copy;
281 std::vector<HGCalGeomParameters::cellParameters> wafers;
283 DDValue val1(attribute, sdTag2, 0.0);
291 <<
" not found but needed.";
293 <<
" not found but needed.";
296 std::unordered_set<std::string>
names;
301 int nsiz = (int)(copy.size());
302 int wafer = (nsiz > 0) ? copy[nsiz-1] : 0;
303 int layer = (nsiz > 1) ? copy[nsiz-2] : 0;
305 edm::LogError(
"HGCalGeom") <<
"Funny wafer # " << wafer <<
" in "
306 << nsiz <<
" components";
309 std::unordered_map<int32_t,int32_t>::iterator itr = copies.find(wafer);
310 std::unordered_map<int32_t,int32_t>::iterator cpy =
311 copiesInLayers[layer].find(wafer);
312 if (itr != copies.end() && cpy == copiesInLayers[layer].end()) {
313 copiesInLayers[layer][wafer] = itr->second;
315 if (itr == copies.end()) {
316 copies[wafer] = wafer2copy.size();
317 copiesInLayers[layer][wafer] = wafer2copy.size();
322 wafer2copy.emplace_back(wafer);
325 wafers.emplace_back(cell);
326 if ( names.count(name) == 0 ) {
328 std::vector<double> zv = polyhedra.
zVec();
329 std::vector<double> rv = polyhedra.
rMaxVec();
331 double dz = 0.5*(zv[1]-zv[0]);
335 mytr.
dz = dz; mytr.
alpha = 0.0;
347 std::map<int,int> wafertype;
348 std::map<int,HGCalGeomParameters::cellParameters> cellsf, cellsc;
349 DDValue val2(attribute, sdTag3, 0.0);
357 <<
" not found but needed.";
359 <<
" not found but needed.";
366 int nsiz = (int)(copy.size());
367 int cellx= (nsiz > 0) ? copy[nsiz-1] : 0;
368 int wafer= (nsiz > 1) ? copy[nsiz-2] : 0;
369 int cell = cellx%1000;
370 int type = cellx/1000;
371 if (type != 1 && type != 2) {
372 edm::LogError(
"HGCalGeom") <<
"Funny cell # " << cell <<
" type "
373 << type <<
" in " << nsiz <<
" components";
376 std::map<int,int>::iterator ktr = wafertype.find(wafer);
377 if (ktr == wafertype.end()) wafertype[wafer] = type;
379 std::map<int,HGCalGeomParameters::cellParameters>::iterator itr;
382 itr = cellsf.find(cell);
383 newc= (itr == cellsf.end());
385 itr = cellsc.find(cell);
386 newc= (itr == cellsc.end());
390 bool half = (name.find(
"Half") != std::string::npos);
399 std::cout <<
"Type " << type <<
" Cell " << cellx <<
" local "
400 << xx <<
":" << yy <<
" new " << p1 <<
":" << p2 <<
"\n";
415 if (((cellsf.size()+cellsc.size())==0) || (wafers.size()==0) ||
416 (layers.size()==0)) {
417 edm::LogError(
"HGCalGeom") <<
"HGCalGeomParameters : number of cells "
418 << cellsf.size() <<
":" << cellsc.size()
419 <<
" wafers " << wafers.size() <<
" layers "
420 << layers.size() <<
" illegal";
422 <<
"HGCalGeomParameters: mismatch between geometry and specpar: cells "
423 << cellsf.size() <<
":" << cellsc.size() <<
" wafers " << wafers.size()
424 <<
" layers " << layers.size();
427 for (
unsigned int i=0;
i<layers.size(); ++
i) {
428 for (std::map<int,HGCalGeomParameters::layerParameters>::iterator itr = layers.begin();
429 itr != layers.end(); ++itr) {
430 if (itr->first == (
int)(
i+1)) {
439 for (
unsigned int i=0;
i<php.
layer_.size(); ++
i) {
440 for (
unsigned int i1=0; i1<trforms.size(); ++i1) {
441 if (!trformUse[i1] && php.
layerGroup_[trforms[i1].lay-1] ==
444 trforms[i1].lay = (
i+1);
445 trformUse[i1] =
true;
448 for (
unsigned int i2=i1+1; i2<trforms.size(); ++i2) {
449 if (!trformUse[i2] && trforms[i2].zp == trforms[i1].zp &&
453 trformUse[i2] =
true;
464 for (
unsigned i = 0;
i < wafer2copy.size(); ++
i ) {
468 std::map<int,int>::iterator ktr = wafertype.find(wafer2copy[i]);
469 int typet = (ktr == wafertype.end()) ? 0 : (ktr->second);
471 double r = wafers[
i].xyz.perp();
473 for (
int k=1;
k<4; ++
k) {
483 std::vector<HGCalGeomParameters::cellParameters>::const_iterator itrf = wafers.end();
484 for (
unsigned int i=0;
i<cellsf.size(); ++
i) {
485 std::map<int,HGCalGeomParameters::cellParameters>::iterator itr = cellsf.find(
i);
486 if (itr == cellsf.end()) {
487 edm::LogError(
"HGCalGeom") <<
"HGCalGeomParameters: missing info for"
488 <<
" fine cell number " <<
i;
490 <<
"HGCalGeomParameters: missing info for fine cell number " <<
i;
492 double xx = (itr->second).xyz.x();
493 double yy = (itr->second).xyz.y();
494 int waf= (itr->second).wafer;
495 std::pair<double,double>
xy =
cellPosition(wafers,itrf,waf,xx,yy);
502 for (
unsigned int i=0;
i<cellsc.size(); ++
i) {
503 std::map<int,HGCalGeomParameters::cellParameters>::iterator itr = cellsc.find(
i);
504 if (itr == cellsc.end()) {
505 edm::LogError(
"HGCalGeom") <<
"HGCalGeomParameters: missing info for"
506 <<
" coarse cell number " <<
i;
508 <<
"HGCalGeomParameters: missing info for coarse cell number " <<
i;
510 double xx = (itr->second).xyz.x();
511 double yy = (itr->second).xyz.y();
512 int waf= (itr->second).wafer;
513 std::pair<double,double>
xy =
cellPosition(wafers,itrf,waf,xx,yy);
547 <<
" layers" << std::endl;
556 <<
" depths" <<std::endl;
559 std::cout <<
"Reco Layer[" <<
i <<
":" << k <<
"] First Layer "
568 <<
" and dimensions of the wafers:" << std::endl;
578 <<
" fine cells in a wafer" << std::endl;
583 <<
" coarse cells in a wafer" << std::endl;
588 <<
" transformation matrices" << std::endl;
603 std::cout <<
"Layer " <<
k <<
":" << theModules.size() << std::endl;
605 for (std::unordered_map<int, int>::const_iterator itr=theModules.begin();
606 itr != theModules.end(); ++itr) {
607 std::cout <<
" " << itr->first <<
":" << itr->second;
609 if (k2 > 9) {
std::cout << std::endl; k2 = 0; }
625 <<
" entries for cellSize_" << std::endl;
634 std::vector<double> slp =
getDDDArray(
"Slope",sv,nmin);
638 <<
" and " << php.
nCells_ <<
" entries for cellFactor_\n";
648 <<
" entries for layerGroup_" << std::endl;
664 for (
unsigned int k=0;
k<php.
boundR_.size(); ++
k)
667 std::cout <<
"HGCalGeomParameters: wafer radius ranges for cell grouping "
673 for (
unsigned int k=0;
k<php.
rLimit_.size(); ++
k)
676 std::cout <<
"HGCalGeomParameters: Minimum/maximum R "
692 std::vector<double> slp =
getDDDArray(
"Slope",sv,nmin);
696 <<
" and layer groupings for the 3 ranges:" << std::endl;
697 for (
int k=0;
k<nmin; ++
k)
704 DDValue val1(attribute, sdTag1, 0.0);
712 std::vector<double> dummy =
getDDDArray(
"WaferSize",sv,nmin);
720 DDValue val2(attribute, sdTag2, 0.0);
732 <<
" cells of sizes:\n";
744 std::cout <<
"Input waferWidth " << waferW <<
":" << rmin
745 <<
" R Limits: " << rin <<
":" << rout
746 <<
" Fine " << rMaxFine <<
"\n";
754 double dx = 0.5*waferW;
755 double dy = 3.0*dx*
tan(30.0*CLHEP::deg);
756 double rr = 2.0*dx*
tan(30.0*CLHEP::deg);
757 int ncol = (int)(2.0*rout/waferW) + 1;
758 int nrow = (int)(rout/(waferW*
tan(30.0*CLHEP::deg))) + 1;
759 int incm(0), inrm(0), kount(0), ntot(0);
763 std::cout <<
"Row " << nrow <<
" Column " << ncol << std::endl;
765 for (
int nr=-nrow; nr <= nrow; ++nr) {
766 int inr = (nr >= 0) ? nr : -nr;
767 for (
int nc=-ncol; nc <= ncol; ++nc) {
768 int inc = (nc >= 0) ? nc : -nc;
769 if (inr%2 == inc%2) {
772 xc[0] = xpos+dx; yc[0] = ypos-0.5*
rr;
773 xc[1] = xpos+dx; yc[1] = ypos+0.5*
rr;
774 xc[2] = xpos; yc[2] = ypos+
rr;
775 xc[3] = xpos-dx; yc[3] = ypos+0.5*
rr;
776 xc[4] = xpos+dx; yc[4] = ypos-0.5*
rr;
777 xc[5] = xpos; yc[5] = ypos-
rr;
778 bool cornerOne(
false);
779 bool cornerAll(
true);
780 for (
int k=0;
k<6; ++
k) {
782 if (rpos >= rin && rpos <= rout) cornerOne =
true;
783 else cornerAll =
false;
785 double rpos =
std::sqrt(xpos*xpos+ypos*ypos);
786 int typet = (rpos < rMaxFine) ? 1 : 2;
788 for (
int k=1;
k<4; ++
k) {
795 int copy = inr*100 + inc;
796 if (nc < 0) copy += 10000;
797 if (nr < 0) copy += 100000;
798 if (inc > incm) incm = inc;
799 if (inr > inrm) inrm = inr;
802 std::cout << kount <<
":" << ntot <<
" Copy " << copy
803 <<
" Type " << typel <<
":" << typet
804 <<
" Location " << cornerOne <<
":" << cornerAll
805 <<
" Position " << xpos <<
":" << ypos <<
"\n";
812 for (
unsigned int il=0; il<php.
layer_.size(); ++il) {
815 for (
int k=0;
k<6; ++
k) {
819 else cornerAll =
false;
822 std::unordered_map<int32_t,int32_t>::iterator cpy =
823 copiesInLayers[php.
layer_[il]].find(copy);
824 if (cpy == copiesInLayers[php.
layer_[il]].end())
825 copiesInLayers[php.
layer_[il]][copy] = cornerAll ? php.
waferCopy_.size() : -1;
835 std::cout <<
"HGCalWaferHexagon: # of columns " << incm <<
" # of rows "
836 << inrm <<
" and " << kount <<
":" << ntot <<
" wafers; R "
837 << rin <<
":" << rout << std::endl;
840 for (
unsigned int k=0;
k<copiesInLayers.size(); ++
k) {
841 const auto& theModules = copiesInLayers[
k];
842 std::cout <<
"Layer " <<
k <<
":" << theModules.size() << std::endl;
844 for (std::unordered_map<int, int>::const_iterator itr=theModules.begin();
845 itr != theModules.end(); ++itr) {
846 std::cout <<
" " << itr->first <<
":" << itr->second;
848 if (k2 > 9) {
std::cout << std::endl; k2 = 0; }
861 DDValue val1(attribute, value, 0.0);
878 <<
" rows for fine cells\n";
883 <<
" rows for coarse cells\n";
895 const std::vector<double> & fvec = value.
doubles();
896 int nval = fvec.size();
899 edm::LogError(
"HGCalGeom") <<
"HGCalGeomParameters : # of " << str
900 <<
" bins " << nval <<
" < " << nmin
902 throw cms::Exception(
"DDException") <<
"HGCalGeomParameters: cannot get array " << str;
905 if (nval < 1 && nmin == 0) {
906 edm::LogError(
"HGCalGeom") <<
"HGCalGeomParameters : # of " << str
907 <<
" bins " << nval <<
" < 1 ==> illegal"
908 <<
" (nmin=" << nmin <<
")";
909 throw cms::Exception(
"DDException") <<
"HGCalGeomParameters: cannot get array " << str;
916 edm::LogError(
"HGCalGeom") <<
"HGCalGeomParameters: cannot get array "
918 throw cms::Exception(
"DDException") <<
"HGCalGeomParameters: cannot get array " << str;
920 std::vector<double> fvec;
926 std::pair<double,double>
928 std::vector<HGCalGeomParameters::cellParameters>::const_iterator& itrf,
929 int wafer,
double xx,
double yy) {
931 if (itrf == wafers.end()) {
932 for (std::vector<HGCalGeomParameters::cellParameters>::const_iterator itr = wafers.begin();
933 itr != wafers.end(); ++itr) {
934 if (itr->wafer == wafer) {
941 if (itrf != wafers.end()) {
942 dx = (xx - itrf->xyz.x());
944 dy = (yy - itrf->xyz.y());
947 return std::pair<double,double>(dx,dy);
std::vector< double > waferPosY_
std::vector< int > layer_
std::vector< double > moduleDzR_
std::vector< int > depthLayerF_
double halfZ(void) const
half of the z-Axis
std::vector< int > depth_
std::vector< double > moduleCellR_
std::vector< double > moduleHR_
const DDLogicalPart & logicalPart() const
The logical-part of the current node in the filtered-view.
const std::vector< double > & doubles() const
a reference to the double-valued values stored in the given instance of DDValue
layer_map copiesInLayers_
double x1(void) const
Half-length along x of the side at y=-pDy1 of the face at -pDz.
std::vector< LayerSetAndLayers > layers(const SeedingLayerSetsHits &sets)
std::vector< bool > cellCoarseHalf_
static const HistoName names[]
std::vector< bool > cellFineHalf_
std::vector< double > rMaxVec(void) const
std::vector< int > moduleLayR_
void addFilter(const DDFilter &, DDLogOp op=DDLogOp::AND)
nav_type copyNumbers() const
return the stack of copy numbers
void loadSpecParsHexagon(const DDFilteredView &, HGCalParameters &, const DDCompactView *, const std::string &, const std::string &)
std::vector< int > cellFine_
std::vector< double > moduleHS_
const DDRotationMatrix & rotation() const
The absolute rotation of the current node.
std::vector< double > trformTranY_
Global3DPoint GlobalPoint
std::vector< double > cellFineY_
std::vector< double > trformRotZY_
const DDSolid & solid(void) const
Returns a reference object of the solid being the shape of this LogicalPart.
std::vector< uint32_t > trformIndex_
void find(edm::Handle< EcalRecHitCollection > &hits, DetId thisDet, std::vector< EcalRecHitCollection::const_iterator > &hit, bool debug=false)
std::vector< int > layerGroupM_
type of data representation of DDCompactView
bool DDfetch(const DDsvalues_type *, DDValue &)
helper for retrieving DDValues from DDsvalues_type *.
std::vector< int > cellFactor_
std::vector< double > trformRotXX_
A DDSolid represents the shape of a part.
void fillTrForm(const hgtrform &mytr)
std::vector< double > trformRotZX_
std::vector< int > dbl_to_int(const std::vector< double > &vecdbl)
Converts a std::vector of doubles to a std::vector of int.
std::vector< double > cellCoarseX_
std::vector< int > cellCoarse_
std::vector< double > trformRotYZ_
std::pair< double, double > cellPosition(const std::vector< cellParameters > &wafers, std::vector< cellParameters >::const_iterator &itrf, int wafer, double xx, double yy)
std::vector< double > boundR_
std::vector< double > cellSize_
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > DD3Vector
A DD Translation is currently implemented with Root Vector3D.
std::vector< double > moduleDzS_
bool next()
set current node to the next node in the filtered tree
std::vector< int > layerIndex_
std::vector< double > moduleAlphaR_
void loadCellParsHexagon(const DDCompactView *cpv, HGCalParameters &php)
std::vector< double > trformRotXY_
Cos< T >::type cos(const T &t)
std::vector< std::pair< unsigned int, DDValue > > DDsvalues_type
std::maps an index to a DDValue. The index corresponds to the index assigned to the name of the std::...
std::vector< double > getDDDArray(const std::string &, const DDsvalues_type &, int &)
std::vector< double > trformRotYX_
hgtrap getModule(unsigned int k, bool reco) const
Interface to a Trapezoid.
Tan< T >::type tan(const T &t)
double y1(void) const
Half-length along y of the face at -pDz.
Abs< T >::type abs(const T &t)
std::vector< std::unordered_map< int32_t, int32_t > > layer_map
std::vector< double > moduleBlR_
std::vector< double > rMinLayHex_
void fillModule(const hgtrap &mytr, bool reco)
std::vector< double > moduleTlS_
std::vector< double > zLayerHex_
std::vector< double > get(const std::string &)
void loadWaferHexagon(HGCalParameters &php)
std::vector< double > rMaxLayHex_
std::vector< double > trformTranX_
std::vector< double > trformRotXZ_
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< double > > XYZPointD
point in space with cartesian internal representation
std::vector< double > zVec(void) const
std::vector< int > layerGroup_
std::vector< double > moduleCellS_
void loadSpecParsSquare(const DDFilteredView &, HGCalParameters &)
double alpha1(void) const
Angle with respect to the y axis from the centre of the side at y=-pDy1 to the centre at y=+pDy1 of t...
DDsvalues_type mergedSpecifics() const
std::vector< double > trformRotYY_
std::vector< double > cellFineX_
double x2(void) const
Half-length along x of the side at y=+pDy1 of the face at -pDz.
std::vector< double > trformRotZZ_
std::vector< double > moduleAlphaS_
std::vector< int > layerGroupO_
std::vector< double > moduleBlS_
bool firstChild()
set the current node to the first child ...
double y2(void) const
Half-length along y of the face at +pDz.
std::vector< int > waferCopy_
std::vector< int > depthIndex_
void loadGeometryHexagon(const DDFilteredView &, HGCalParameters &, const std::string &, const DDCompactView *, const std::string &, const std::string &)
std::vector< double > rLimit_
std::vector< int > waferTypeT_
const DDTranslation & translation() const
The absolute translation of the current node.
std::vector< double > cellCoarseY_
std::vector< int > moduleLayS_
void loadGeometrySquare(const DDFilteredView &, HGCalParameters &, const std::string &)
std::vector< double > trformTranZ_
for(const auto &isodef:isoDefs)
std::vector< double > waferPosX_
void addTrForm(const CLHEP::Hep3Vector &h3v)
std::vector< double > moduleTlR_
std::vector< int > waferTypeL_
void setCriteria(const DDValue &nameVal, DDCompOp, DDLogOp l=DDLogOp::AND, bool asString=true, bool merged=true)
The DDGenericFilter is a runtime-parametrized Filter looking on DDSpecifcs.