Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:14:51 

0001 #include "Geometry/HcalTowerAlgo/interface/HcalFlexiHardcodeGeometryLoader.h"
0002 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0003 #include "Geometry/CaloGeometry/interface/IdealObliquePrism.h"
0004 #include "Geometry/CaloGeometry/interface/IdealZPrism.h"
0005 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0006 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0007 #include "Geometry/HcalCommonData/interface/HcalDDDRecConstants.h"
0008 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0009 
0010 #include <algorithm>
0011 #include <sstream>
0012 #include <vector>
0013 
0014 using CCGFloat = CaloCellGeometry::CCGFloat;
0015 
0016 //#define EDM_ML_DEBUG
0017 // ==============> Loader Itself <==========================
0018 
0019 HcalFlexiHardcodeGeometryLoader::HcalFlexiHardcodeGeometryLoader() {
0020   MAX_HCAL_PHI = 72;
0021   DEGREE2RAD = M_PI / 180.;
0022 }
0023 
0024 CaloSubdetectorGeometry* HcalFlexiHardcodeGeometryLoader::load(const HcalTopology& fTopology,
0025                                                                const HcalDDDRecConstants& hcons) {
0026   HcalGeometry* hcalGeometry = new HcalGeometry(fTopology);
0027   if (nullptr == hcalGeometry->cornersMgr())
0028     hcalGeometry->allocateCorners(fTopology.ncells() + fTopology.getHFSize());
0029   if (nullptr == hcalGeometry->parMgr())
0030     hcalGeometry->allocatePar(hcalGeometry->numberOfShapes(), HcalGeometry::k_NumberOfParametersPerShape);
0031   isBH_ = hcons.isBH();
0032 #ifdef EDM_ML_DEBUG
0033   edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader initialize with ncells " << fTopology.ncells() << " and shapes "
0034                                << hcalGeometry->numberOfShapes() << ":" << HcalGeometry::k_NumberOfParametersPerShape
0035                                << " with BH Flag " << isBH_;
0036 #endif
0037   if (fTopology.mode() == HcalTopologyMode::H2) {  // TB geometry
0038     fillHBHO(hcalGeometry, makeHBCells(hcons), true);
0039     fillHBHO(hcalGeometry, makeHOCells(), false);
0040     fillHE(hcalGeometry, makeHECells_H2());
0041   } else {  // regular geometry
0042     fillHBHO(hcalGeometry, makeHBCells(hcons), true);
0043     fillHBHO(hcalGeometry, makeHOCells(), false);
0044     fillHF(hcalGeometry, makeHFCells(hcons));
0045     fillHE(hcalGeometry, makeHECells(hcons));
0046   }
0047   //fast insertion of valid ids requires sort at end
0048   hcalGeometry->sortValidIds();
0049   return hcalGeometry;
0050 }
0051 
0052 // ----------> HB <-----------
0053 std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> HcalFlexiHardcodeGeometryLoader::makeHBCells(
0054     const HcalDDDRecConstants& hcons) {
0055   std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> result;
0056   std::vector<std::pair<double, double> > gconsHB = hcons.getConstHBHE(0);
0057   std::vector<HcalDDDRecConstants::HcalEtaBin> etabins = hcons.getEtaBins(0);
0058 
0059 #ifdef EDM_ML_DEBUG
0060   edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader called for " << etabins.size() << " Eta Bins";
0061   for (unsigned int k = 0; k < gconsHB.size(); ++k) {
0062     edm::LogVerbatim("HCalGeom") << "gconsHB[" << k << "] = " << gconsHB[k].first << " +- " << gconsHB[k].second;
0063   }
0064 #endif
0065   for (auto& etabin : etabins) {
0066     int iring = (etabin.zside >= 0) ? etabin.ieta : -etabin.ieta;
0067     int depth = etabin.depthStart;
0068     double dphi =
0069         (etabin.phis.size() > 1) ? (etabin.phis[1].second - etabin.phis[0].second) : ((2.0 * M_PI) / MAX_HCAL_PHI);
0070     for (unsigned int k = 0; k < etabin.layer.size(); ++k) {
0071       int layf = etabin.layer[k].first - 1;
0072       int layl = etabin.layer[k].second - 1;
0073       double rmin = gconsHB[layf].first - gconsHB[layf].second;
0074       double rmax = gconsHB[layl].first + gconsHB[layl].second;
0075       for (unsigned int j = 0; j < etabin.phis.size(); ++j) {
0076 #ifdef EDM_ML_DEBUG
0077         edm::LogVerbatim("HCalGeom") << "HBRing " << iring << " eta " << etabin.etaMin << ":" << etabin.etaMax
0078                                      << " depth " << depth << " R " << rmin << ":" << rmax << " Phi "
0079                                      << etabin.phis[j].first << ":" << etabin.phis[j].second << ":" << dphi << " layer["
0080                                      << k << "]: " << etabin.layer[k].first - 1 << ":" << etabin.layer[k].second;
0081 #endif
0082         result.emplace_back(HcalFlexiHardcodeGeometryLoader::HBHOCellParameters(
0083             iring, depth, etabin.phis[j].first, etabin.phis[j].second, dphi, rmin, rmax, etabin.etaMin, etabin.etaMax));
0084       }
0085       depth++;
0086     }
0087   }
0088   return result;
0089 }
0090 
0091 // ----------> HO <-----------
0092 std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> HcalFlexiHardcodeGeometryLoader::makeHOCells() {
0093   const double HORMIN0 = 390.0;
0094   const double HORMIN1 = 412.6;
0095   const double HORMAX = 413.6;
0096   const int nCells = 15;
0097   const int nPhi = 72;
0098   const double etamin[nCells] = {
0099       0.000, 0.087, 0.174, 0.261, 0.3395, 0.435, 0.522, 0.609, 0.696, 0.783, 0.873, 0.957, 1.044, 1.131, 1.218};
0100   const double etamax[nCells] = {
0101       0.087, 0.174, 0.261, 0.3075, 0.435, 0.522, 0.609, 0.696, 0.783, 0.8494, 0.957, 1.044, 1.131, 1.218, 1.305};
0102   std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters> result;
0103   result.reserve(nCells * nPhi);
0104   double dphi = ((2.0 * M_PI) / nPhi);
0105   for (int i = 0; i < nCells; ++i) {
0106     double rmin = ((i < 4) ? HORMIN0 : HORMIN1);
0107     for (int iside = -1; iside <= 1; iside += 2) {
0108       for (int j = 0; j < nPhi; ++j) {
0109         double phi = (j + 0.5) * dphi;
0110         // eta, depth, phi, phi0, deltaPhi, rMin, rMax, etaMin, etaMax
0111         result.emplace_back(HcalFlexiHardcodeGeometryLoader::HBHOCellParameters(
0112             iside * (i + 1), 4, j + 1, phi, dphi, rmin, HORMAX, etamin[i], etamax[i]));
0113       }
0114     }
0115   }
0116   return result;
0117 }
0118 
0119 //
0120 // Convert constants to appropriate cells
0121 //
0122 void HcalFlexiHardcodeGeometryLoader::fillHBHO(
0123     HcalGeometry* fGeometry, const std::vector<HcalFlexiHardcodeGeometryLoader::HBHOCellParameters>& fCells, bool fHB) {
0124   fGeometry->increaseReserve(fCells.size());
0125   for (const auto& param : fCells) {
0126     HcalDetId hid(fHB ? HcalBarrel : HcalOuter, param.ieta, param.iphi, param.depth);
0127     float phiCenter = param.phi;  // middle of the cell
0128     float etaCenter = 0.5 * (param.etaMin + param.etaMax);
0129     float x = param.rMin * cos(phiCenter);
0130     float y = param.rMin * sin(phiCenter);
0131     float z = (param.ieta < 0) ? -(param.rMin * sinh(etaCenter)) : (param.rMin * sinh(etaCenter));
0132     // make cell geometry
0133     GlobalPoint refPoint(x, y, z);  // center of the cell's face
0134     std::vector<CCGFloat> cellParams;
0135     cellParams.reserve(5);
0136     cellParams.emplace_back(0.5 * (param.etaMax - param.etaMin));                // deta_half
0137     cellParams.emplace_back(0.5 * param.dphi);                                   // dphi_half
0138     cellParams.emplace_back(0.5 * (param.rMax - param.rMin) * cosh(etaCenter));  // dr_half
0139     cellParams.emplace_back(fabs(refPoint.eta()));
0140     cellParams.emplace_back(fabs(refPoint.z()));
0141 #ifdef EDM_ML_DEBUG
0142     edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader::fillHBHO-> " << hid << " " << hid.rawId() << " "
0143                                  << std::hex << hid.rawId() << std::dec << " " << hid << " " << refPoint << '/'
0144                                  << cellParams[0] << '/' << cellParams[1] << '/' << cellParams[2];
0145 #endif
0146     fGeometry->newCellFast(refPoint,
0147                            refPoint,
0148                            refPoint,
0149                            CaloCellGeometry::getParmPtr(cellParams, fGeometry->parMgr(), fGeometry->parVecVec()),
0150                            hid);
0151   }
0152 }
0153 
0154 // ----------> HE <-----------
0155 std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> HcalFlexiHardcodeGeometryLoader::makeHECells(
0156     const HcalDDDRecConstants& hcons) {
0157   std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> result;
0158   std::vector<std::pair<double, double> > gconsHE = hcons.getConstHBHE(1);
0159 #ifdef EDM_ML_DEBUG
0160   edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader:HE with " << gconsHE.size() << " cells";
0161 #endif
0162   if (!gconsHE.empty()) {
0163     std::vector<HcalDDDRecConstants::HcalEtaBin> etabins = hcons.getEtaBins(1);
0164 
0165 #ifdef EDM_ML_DEBUG
0166     edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader called for HE with " << etabins.size() << " Eta Bins and "
0167                                  << gconsHE.size() << " depths";
0168     std::ostringstream st1;
0169     for (unsigned int i = 0; i < gconsHE.size(); ++i)
0170       st1 << " Depth[" << i << "] = " << gconsHE[i].first << " +- " << gconsHE[i].second;
0171     edm::LogVerbatim("HCalGeom") << st1.str();
0172 #endif
0173     for (auto& etabin : etabins) {
0174       int iring = (etabin.zside >= 0) ? etabin.ieta : -etabin.ieta;
0175       int depth = etabin.depthStart;
0176       double dphi =
0177           (etabin.phis.size() > 1) ? (etabin.phis[1].second - etabin.phis[0].second) : ((2.0 * M_PI) / MAX_HCAL_PHI);
0178 #ifdef EDM_ML_DEBUG
0179       edm::LogVerbatim("HCalGeom") << "FlexiGeometryLoader::Ring " << iring << " nphi " << etabin.phis.size()
0180                                    << " dstart " << depth << " dphi " << dphi << " layers " << etabin.layer.size();
0181       std::ostringstream st2;
0182       for (unsigned int j = 0; j < etabin.phis.size(); ++j)
0183         st2 << " [" << j << "] " << etabin.phis[j].first << ":" << etabin.phis[j].second;
0184       edm::LogVerbatim("HCalGeom") << st2.str();
0185 #endif
0186       for (unsigned int k = 0; k < etabin.layer.size(); ++k) {
0187         int layf = etabin.layer[k].first - 1;
0188         int layl = etabin.layer[k].second - 1;
0189         double zmin = gconsHE[layf].first - gconsHE[layf].second;
0190         double zmax = gconsHE[layl].first + gconsHE[layl].second;
0191         if (zmin < 1.0) {
0192           for (int k2 = layf; k2 <= layl; ++k2) {
0193             if (gconsHE[k2].first > 10) {
0194               zmin = gconsHE[k2].first - gconsHE[k2].second;
0195               break;
0196             }
0197           }
0198         }
0199         if (zmin >= zmax)
0200           zmax = zmin + 10.;
0201         for (unsigned int j = 0; j < etabin.phis.size(); ++j) {
0202 #ifdef EDM_ML_DEBUG
0203           edm::LogVerbatim("HCalGeom") << "HERing " << iring << " eta " << etabin.etaMin << ":" << etabin.etaMax
0204                                        << " depth " << depth << " Z " << zmin << ":" << zmax
0205                                        << " Phi :" << etabin.phis[j].first << ":" << etabin.phis[j].second << ":"
0206                                        << dphi << " layer[" << k << "]: " << etabin.layer[k].first - 1 << ":"
0207                                        << etabin.layer[k].second - 1;
0208 #endif
0209           result.emplace_back(HcalFlexiHardcodeGeometryLoader::HECellParameters(
0210               iring, depth, etabin.phis[j].first, etabin.phis[j].second, dphi, zmin, zmax, etabin.etaMin, etabin.etaMax));
0211         }
0212         depth++;
0213       }
0214     }
0215   }
0216   return result;
0217 }
0218 
0219 // ----------> HE @ H2 <-----------
0220 std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> HcalFlexiHardcodeGeometryLoader::makeHECells_H2() {
0221   const double HEZMIN_H2 = 400.715;
0222   const double HEZMID_H2 = 436.285;
0223   const double HEZMAX_H2 = 541.885;
0224   const int nEtas = 10;
0225   const int nDepth[nEtas] = {1, 2, 2, 2, 2, 2, 2, 2, 3, 3};
0226   const int dStart[nEtas] = {3, 1, 1, 1, 1, 1, 1, 1, 1, 1};
0227   const int nPhis[nEtas] = {8, 8, 8, 8, 8, 8, 4, 4, 4, 4};
0228   const double etas[nEtas + 1] = {1.305, 1.373, 1.444, 1.521, 1.603, 1.693, 1.790, 1.880, 1.980, 2.090, 2.210};
0229   const double zval[4 * nEtas] = {
0230       409.885,   462.685,   0.,        0.,        HEZMIN_H2, 427.485,   506.685,   0.0,       HEZMIN_H2, HEZMID_H2,
0231       524.285,   0.,        HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0.,        HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0.,
0232       HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0.,        HEZMIN_H2, HEZMID_H2, HEZMAX_H2, 0.,        HEZMIN_H2, HEZMID_H2,
0233       HEZMAX_H2, 0.,        HEZMIN_H2, 418.685,   HEZMID_H2, HEZMAX_H2, HEZMIN_H2, 418.685,   HEZMID_H2, HEZMAX_H2};
0234   std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters> result;
0235 
0236   for (int i = 0; i < nEtas; ++i) {
0237     int ieta = 16 + i;
0238     for (int k = 0; k < nDepth[i]; ++k) {
0239       int depth = dStart[i] + k;
0240       for (int j = 0; j < nPhis[i]; ++j) {
0241         int iphi = (nPhis[i] == 8) ? (j + 1) : (2 * j + 1);
0242         double dphi = (40.0 * DEGREE2RAD) / nPhis[i];
0243         double phi0 = (j + 0.5) * dphi;
0244         // ieta, depth, iphi, phi0, deltaPhi, zMin, zMax, etaMin, etaMax
0245         result.emplace_back(HcalFlexiHardcodeGeometryLoader::HECellParameters(
0246             ieta, depth, iphi, phi0, dphi, zval[4 * i + k + 1], zval[4 * i + k + 2], etas[i], etas[i + 1]));
0247       }
0248     }
0249   }
0250   return result;
0251 }
0252 
0253 // ----------> HF <-----------
0254 std::vector<HcalFlexiHardcodeGeometryLoader::HFCellParameters> HcalFlexiHardcodeGeometryLoader::makeHFCells(
0255     const HcalDDDRecConstants& hcons) {
0256   const float HFZMIN1 = 1115.;
0257   const float HFZMIN2 = 1137.;
0258   const float HFZMAX = 1280.1;
0259   std::vector<HcalDDDRecConstants::HFCellParameters> cells = hcons.getHFCellParameters();
0260   unsigned int nCells = cells.size();
0261   std::vector<HcalFlexiHardcodeGeometryLoader::HFCellParameters> result;
0262   result.reserve(nCells);
0263   for (unsigned int i = 0; i < nCells; ++i) {
0264     HcalFlexiHardcodeGeometryLoader::HFCellParameters cell1(cells[i].ieta,
0265                                                             cells[i].depth,
0266                                                             cells[i].firstPhi,
0267                                                             cells[i].stepPhi,
0268                                                             cells[i].nPhi,
0269                                                             5 * cells[i].stepPhi,
0270                                                             HFZMIN1,
0271                                                             HFZMAX,
0272                                                             cells[i].rMin,
0273                                                             cells[i].rMax);
0274     result.emplace_back(cell1);
0275     HcalFlexiHardcodeGeometryLoader::HFCellParameters cell2(cells[i].ieta,
0276                                                             1 + cells[i].depth,
0277                                                             cells[i].firstPhi,
0278                                                             cells[i].stepPhi,
0279                                                             cells[i].nPhi,
0280                                                             5 * cells[i].stepPhi,
0281                                                             HFZMIN2,
0282                                                             HFZMAX,
0283                                                             cells[i].rMin,
0284                                                             cells[i].rMax);
0285     result.emplace_back(cell2);
0286   }
0287   return result;
0288 }
0289 
0290 void HcalFlexiHardcodeGeometryLoader::fillHE(
0291     HcalGeometry* fGeometry, const std::vector<HcalFlexiHardcodeGeometryLoader::HECellParameters>& fCells) {
0292   fGeometry->increaseReserve(fCells.size());
0293   for (const auto& param : fCells) {
0294     HcalDetId hid(HcalEndcap, param.ieta, param.iphi, param.depth);
0295     float phiCenter = param.phi;  // middle of the cell
0296     float etaCenter = 0.5 * (param.etaMin + param.etaMax);
0297     int iside = (param.ieta >= 0) ? 1 : -1;
0298     float z = (isBH_) ? (iside * 0.5 * (param.zMin + param.zMax)) : (iside * param.zMin);
0299     float perp = std::abs(z) / sinh(etaCenter);
0300     float x = perp * cos(phiCenter);
0301     float y = perp * sin(phiCenter);
0302     // make cell geometry
0303     GlobalPoint refPoint(x, y, z);  // center of the cell's face
0304     std::vector<CCGFloat> cellParams;
0305     cellParams.reserve(5);
0306     cellParams.emplace_back(0.5 * (param.etaMax - param.etaMin));                 //deta_half
0307     cellParams.emplace_back(0.5 * param.dphi);                                    // dphi_half
0308     cellParams.emplace_back(-0.5 * (param.zMax - param.zMin) / tanh(etaCenter));  // dz_half, "-" means edges in Z
0309     cellParams.emplace_back(fabs(refPoint.eta()));
0310     cellParams.emplace_back(fabs(refPoint.z()));
0311 #ifdef EDM_ML_DEBUG
0312     edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader::fillHE-> " << hid << " " << hid.rawId() << " "
0313                                  << std::hex << hid.rawId() << std::dec << " " << hid << refPoint << '/'
0314                                  << cellParams[0] << '/' << cellParams[1] << '/' << cellParams[2];
0315 #endif
0316     fGeometry->newCellFast(refPoint,
0317                            refPoint,
0318                            refPoint,
0319                            CaloCellGeometry::getParmPtr(cellParams, fGeometry->parMgr(), fGeometry->parVecVec()),
0320                            hid);
0321   }
0322 }
0323 
0324 void HcalFlexiHardcodeGeometryLoader::fillHF(
0325     HcalGeometry* fGeometry, const std::vector<HcalFlexiHardcodeGeometryLoader::HFCellParameters>& fCells) {
0326   fGeometry->increaseReserve(fCells.size());
0327   for (const auto& param : fCells) {
0328     for (int kPhi = 0; kPhi < param.nPhi; ++kPhi) {
0329       int iPhi = param.phiFirst + kPhi * param.phiStep;
0330       HcalDetId hid(HcalForward, param.eta, iPhi, param.depth);
0331       // middle of the cell
0332       float phiCenter = ((iPhi - 1) * 360. / MAX_HCAL_PHI + 0.5 * param.dphi) * DEGREE2RAD;
0333       GlobalPoint inner(param.rMin, 0, param.zMin);
0334       GlobalPoint outer(param.rMax, 0, param.zMin);
0335       float iEta = inner.eta();
0336       float oEta = outer.eta();
0337       float etaCenter = 0.5 * (iEta + oEta);
0338 
0339       float perp = param.zMin / sinh(etaCenter);
0340       float x = perp * cos(phiCenter);
0341       float y = perp * sin(phiCenter);
0342       float z = (param.eta > 0) ? param.zMin : -param.zMin;
0343       // make cell geometry
0344       GlobalPoint refPoint(x, y, z);  // center of the cell's face
0345       std::vector<CCGFloat> cellParams;
0346       cellParams.reserve(5);
0347       cellParams.emplace_back(0.5 * (iEta - oEta));              // deta_half
0348       cellParams.emplace_back(0.5 * param.dphi * DEGREE2RAD);    // dphi_half
0349       cellParams.emplace_back(0.5 * (param.zMax - param.zMin));  // dz_half
0350       cellParams.emplace_back(fabs(refPoint.eta()));
0351       cellParams.emplace_back(fabs(refPoint.z()));
0352 #ifdef EDM_ML_DEBUG
0353       edm::LogVerbatim("HCalGeom") << "HcalFlexiHardcodeGeometryLoader::fillHF-> " << hid << " " << hid.rawId() << " "
0354                                    << std::hex << hid.rawId() << std::dec << " " << hid << " " << refPoint << '/'
0355                                    << cellParams[0] << '/' << cellParams[1] << '/' << cellParams[2];
0356 #endif
0357       fGeometry->newCellFast(refPoint,
0358                              refPoint,
0359                              refPoint,
0360                              CaloCellGeometry::getParmPtr(cellParams, fGeometry->parMgr(), fGeometry->parVecVec()),
0361                              hid);
0362     }
0363   }
0364 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!