Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-14 22:36:02

0001 
0002 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0003 
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/Framework/interface/MakerMacros.h"
0006 #include "FWCore/Utilities/interface/ESGetToken.h"
0007 
0008 #include "DataFormats/HcalDetId/interface/HcalCastorDetId.h"
0009 #include "DataFormats/HcalDetId/interface/HcalZDCDetId.h"
0010 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0011 
0012 #include "DataFormats/Math/interface/GeantUnits.h"
0013 
0014 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0015 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
0016 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0017 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0018 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0019 #include "Geometry/CaloGeometry/interface/TruncatedPyramid.h"
0020 #include "Geometry/EcalAlgo/interface/EcalPreshowerGeometry.h"
0021 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0022 #include "Geometry/CaloGeometry/interface/CaloGenericDetId.h"
0023 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0024 #include "Geometry/CaloTopology/interface/CaloSubdetectorTopology.h"
0025 
0026 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0027 
0028 #include "FWCore/ServiceRegistry/interface/Service.h"
0029 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0030 
0031 #include <fstream>
0032 #include <iomanip>
0033 #include <iterator>
0034 #include "TH1.h"
0035 #include "TH1D.h"
0036 #include "TProfile.h"
0037 
0038 using namespace geant_units;
0039 using namespace geant_units::operators;
0040 
0041 class CaloGeometryAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0042   enum CenterOrCorner { kCenter, kCorner };
0043   enum XorYorZ { kX, kY, kZ };
0044 
0045 public:
0046   explicit CaloGeometryAnalyzer(const edm::ParameterSet&);
0047   ~CaloGeometryAnalyzer() override;
0048 
0049   void beginJob() override {}
0050   void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0051   void endJob() override {}
0052 
0053 private:
0054   // ----------member data ---------------------------
0055   void build(const CaloGeometry* cg,
0056              const HcalTopology& ht,
0057              DetId::Detector det,
0058              int subdetn,
0059              const char* name,
0060              unsigned int histi);
0061 
0062   void buildHcal(const CaloGeometry* cg,
0063                  const HcalTopology& ht,
0064                  DetId::Detector det,
0065                  int subdetn,
0066                  const char* name,
0067                  unsigned int histi);
0068 
0069   void ctrcor(const DetId& did,
0070               const CaloCellGeometry& cell,
0071               std::fstream& fCtr,
0072               std::fstream& fCor,
0073               std::fstream& oldCtr,
0074               std::fstream& oldCor,
0075               unsigned int histi);
0076 
0077   void checkDiff(int i1, int i2, int i3, CenterOrCorner iCtrCor, XorYorZ iXYZ, double diff);
0078   int pass_;
0079 
0080   EEDetId gid(unsigned int ix, unsigned int iy, unsigned int iz, const EEDetId& did) const;
0081 
0082   void cmpset(const CaloSubdetectorGeometry* geom, const GlobalPoint& gp, const double dR);
0083 
0084   void ovrTst(const CaloGeometry* cg, const CaloSubdetectorGeometry* geom, const EEDetId& id, std::fstream& fOvr);
0085 
0086   void ovrTst(const CaloGeometry* cg, const CaloSubdetectorGeometry* geom, const EBDetId& id, std::fstream& fOvr);
0087 
0088   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometryToken_;
0089   edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> topologyToken_;
0090 
0091   edm::Service<TFileService> h_fs;
0092 
0093   TProfile* h_dPhi[7];
0094   TProfile* h_dPhiR[7];
0095 
0096   TProfile* h_dEta[7];
0097   TProfile* h_dEtaR[7];
0098 
0099   TProfile* h_eta;
0100   TProfile* h_phi;
0101 
0102   TH1D* h_diffs[10][12];
0103 
0104   TH1D* h_scindex;
0105 
0106   bool m_allOK;
0107 };
0108 
0109 CaloGeometryAnalyzer::CaloGeometryAnalyzer(const edm::ParameterSet& /*iConfig*/)
0110     : geometryToken_{esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{})},
0111       topologyToken_{esConsumes<HcalTopology, HcalRecNumberingRecord>(edm::ESInputTag{})} {
0112   usesResource("TFileService");
0113 
0114   pass_ = 0;
0115 
0116   h_dPhi[0] = h_fs->make<TProfile>("dPhi:EB:index", "EB: dPhi vs index", 61200, -0.5, 61199.5, " ");
0117   h_dPhiR[0] = h_fs->make<TProfile>("dPhi:EB:R", "EB: dPhi vs R", 100, 125, 135, " ");
0118 
0119   h_dEta[0] = h_fs->make<TProfile>("dEta:EB:index", "EB: dEta vs index", 61200, -0.5, 61199.5, " ");
0120   h_dEtaR[0] = h_fs->make<TProfile>("dEta:EB:R", "EB: dEta vs R", 100, 125, 135, " ");
0121 
0122   h_dPhi[1] = h_fs->make<TProfile>("dPhi:EE:index", "EE: dPhi vs index", 14648, -0.5, 14647.5, " ");
0123   h_dPhiR[1] = h_fs->make<TProfile>("dPhi:EE:R", "EE: dPhi vs R", 130, 30, 160, " ");
0124 
0125   h_dEta[1] = h_fs->make<TProfile>("dEta:EE:index", "EE: dEta vs index", 14648, -0.5, 14647.5, " ");
0126   h_dEtaR[1] = h_fs->make<TProfile>("dEta:EE:R", "EE: dEta vs R", 130, 30, 160, " ");
0127 
0128   h_dPhi[2] = h_fs->make<TProfile>("dPhi:ES:index", "ES: dPhi vs index", 137216, -0.5, 137215.5, " ");
0129   h_dPhiR[2] = h_fs->make<TProfile>("dPhi:ES:R", "ES: dPhi vs R", 90, 40, 130, " ");
0130 
0131   h_dEta[2] = h_fs->make<TProfile>("dEta:ES:index", "ES: dEta vs index", 137216, -0.5, 137215.5, " ");
0132   h_dEtaR[2] = h_fs->make<TProfile>("dEta:ES:R", "ES: dEta vs R", 90, 40, 130, " ");
0133 
0134   h_dPhi[3] = h_fs->make<TProfile>("dPhi:HC:index", "HC: dPhi vs index", 9072, -0.5, 9071.5, " ");
0135   h_dPhiR[3] = h_fs->make<TProfile>("dPhi:HC:R", "HC: dPhi vs R", 400, 0, 400, " ");
0136 
0137   h_dEta[3] = h_fs->make<TProfile>("dEta:HC:index", "HC: dEta vs index", 9072, -0.5, 9071.5, " ");
0138   h_dEtaR[3] = h_fs->make<TProfile>("dEta:HC:R", "HC: dEta vs R", 400, 0, 400, " ");
0139 
0140   h_dPhi[4] = h_fs->make<TProfile>("dPhi:ZD:index", "ZD: dPhi vs index", 22, -0.5, 21.5, " ");
0141   h_dPhiR[4] = h_fs->make<TProfile>("dPhi:ZD:R", "ZD: dPhi vs R", 100, 0, 10, " ");
0142 
0143   h_dEta[4] = h_fs->make<TProfile>("dEta:ZD:index", "ZD: dEta vs index", 22, -0.5, 21.5, " ");
0144   h_dEtaR[4] = h_fs->make<TProfile>("dEta:ZD:R", "ZD: dEta vs R", 100, 0, 10, " ");
0145 
0146   h_dPhi[5] = h_fs->make<TProfile>("dPhi:CA:index", "CA: dPhi vs index", 224, -0.5, 223.5, " ");
0147   h_dPhiR[5] = h_fs->make<TProfile>("dPhi:CA:R", "CA: dPhi vs R", 100, 0, 20, " ");
0148 
0149   h_dEta[5] = h_fs->make<TProfile>("dEta:CA:index", "CA: dEta vs index", 224, -0.5, 223.5, " ");
0150   h_dEtaR[5] = h_fs->make<TProfile>("dEta:CA:R", "CA: dEta vs R", 100, 0, 20, " ");
0151 
0152   h_dPhi[6] = h_fs->make<TProfile>("dPhi:CT:index", "CT: dPhi vs index", 4320, -0.5, 4319.5, " ");
0153   h_dPhiR[6] = h_fs->make<TProfile>("dPhi:CT:R", "CT: dPhi vs R", 150, 0, 150, " ");
0154 
0155   h_dEta[6] = h_fs->make<TProfile>("dEta:CT:index", "CT: dEta vs index", 4320, -0.5, 4319.5, " ");
0156   h_dEtaR[6] = h_fs->make<TProfile>("dEta:CT:R", "CT: dEta vs R", 150, 0, 150, " ");
0157 
0158   h_eta = h_fs->make<TProfile>("iEta", "Eta vs iEta", 86 * 2 * 4, -86, 86, " ");
0159   h_phi = h_fs->make<TProfile>("iPhi", "Phi vs iPhi", 360 * 4, 1, 361, " ");
0160 
0161   const std::string hname[10] = {"EB", "EE", "ES", "HB", "HO", "HE", "HF", "CT", "ZD", "CA"};
0162   const std::string cname[12] = {
0163       "XCtr", "YCtr", "ZCtr", "XCor0", "YCor0", "ZCor0", "XCor3", "YCor3", "ZCor3", "XCor6", "YCor6", "ZCor6"};
0164 
0165   for (unsigned int i(0); i != 10; ++i) {
0166     for (unsigned int j(0); j != 12; ++j) {
0167       h_diffs[i][j] =
0168           h_fs->make<TH1D>(std::string(hname[i] + cname[j] + std::string("Diff (microns)")).c_str(),
0169                            std::string(hname[i] + std::string(": New-Nom(") + cname[j] + std::string(")")).c_str(),
0170                            200,
0171                            -200.,
0172                            200.);
0173     }
0174   }
0175   h_scindex = h_fs->make<TH1D>(
0176       std::string("Supercrystal Hashed Index").c_str(), std::string("SC Hashed Index").c_str(), 632, -0.5, 631.5);
0177 }
0178 
0179 CaloGeometryAnalyzer::~CaloGeometryAnalyzer() {}
0180 
0181 void CaloGeometryAnalyzer::cmpset(const CaloSubdetectorGeometry* geom, const GlobalPoint& gp, const double dR) {
0182   typedef CaloSubdetectorGeometry::DetIdSet DetSet;
0183 
0184   DetSet base = geom->CaloSubdetectorGeometry::getCells(gp, dR);
0185   DetSet over = geom->getCells(gp, dR);
0186   if (over == base) {
0187     /*
0188       std::cout << "getCells Test dR="
0189         << dR
0190         << ", gp=" << gp 
0191         << ", gp.eta=" << gp.eta() 
0192         << ", gp.phi=" << gp.phi() 
0193         << ": base and over are equal!\n ***************************\n " 
0194         << std::endl ;
0195 */
0196   } else {
0197     if (2 < std::abs((int)(base.size()) - (int)(over.size()))) {
0198       DetSet inBaseNotOver;
0199       DetSet inOverNotBase;
0200       std::set_difference(
0201           base.begin(), base.end(), over.begin(), over.end(), std::inserter(inBaseNotOver, inBaseNotOver.begin()));
0202 
0203       if (inBaseNotOver.empty()) {
0204         std::cout << "getCells Test dR=" << dR << ", gp=" << gp << ", gp.eta=" << gp.eta() << ", gp.phi=" << gp.phi()
0205                   << ": No elements in base but not overload " << std::endl;
0206       } else {
0207         std::cout << "Length of Base is " << base.size() << std::endl;
0208         std::cout << "Length of Over is " << over.size() << std::endl;
0209         std::cout << "There are " << inBaseNotOver.size() << " items in Base but not in Overload" << std::endl;
0210 
0211         for (const auto& iS : inBaseNotOver) {
0212           std::cout << "getCells Test dR=" << dR << ", gp=" << gp << ": cell in base but not overload = ";
0213           if (iS.det() == DetId::Ecal && iS.subdetId() == EcalBarrel)
0214             std::cout << EBDetId(iS);
0215           if (iS.det() == DetId::Ecal && iS.subdetId() == EcalEndcap)
0216             std::cout << EEDetId(iS);
0217           if (iS.det() == DetId::Hcal)
0218             std::cout << HcalDetId(iS);
0219           std::cout << std::endl;
0220         }
0221       }
0222 
0223       std::set_difference(
0224           over.begin(), over.end(), base.begin(), base.end(), std::inserter(inOverNotBase, inOverNotBase.begin()));
0225 
0226       if (inOverNotBase.empty()) {
0227         std::cout << "getCells Test dR=" << dR << ", gp=" << gp << ", gp.eta=" << gp.eta() << ", gp.phi=" << gp.phi()
0228                   << ": No elements in overload but not base " << std::endl;
0229       } else {
0230         std::cout << "Length of Base is " << base.size() << std::endl;
0231         std::cout << "Length of Over is " << over.size() << std::endl;
0232         std::cout << "There are " << inOverNotBase.size() << " items in Overload but not in Base" << std::endl;
0233 
0234         for (const auto& iS : inOverNotBase) {
0235           std::cout << "getCells Test dR=" << dR << ", gp=" << gp << ": cell in overload but not base = ";
0236           if (iS.det() == DetId::Ecal && iS.subdetId() == EcalBarrel)
0237             std::cout << EBDetId(iS);
0238           if (iS.det() == DetId::Ecal && iS.subdetId() == EcalEndcap)
0239             std::cout << EEDetId(iS);
0240           if (iS.det() == DetId::Hcal)
0241             std::cout << HcalDetId(iS);
0242           std::cout << std::endl;
0243         }
0244       }
0245       std::cout << "------------- done with mismatch printout ---------------" << std::endl;
0246     }
0247   }
0248 }
0249 
0250 EEDetId CaloGeometryAnalyzer::gid(unsigned int ix, unsigned int iy, unsigned int iz, const EEDetId& did) const {
0251   return (EEDetId::validDetId(ix, iy, iz)
0252               ? EEDetId(ix, iy, iz)
0253               : (EEDetId::validDetId(ix + 1, iy, iz)
0254                      ? EEDetId(ix + 1, iy, iz)
0255                      : (EEDetId::validDetId(ix - 1, iy, iz)
0256                             ? EEDetId(ix - 1, iy, iz)
0257                             : (EEDetId::validDetId(ix, iy + 1, iz)
0258                                    ? EEDetId(ix, iy + 1, iz)
0259                                    : (EEDetId::validDetId(ix, iy - 1, iz)
0260                                           ? EEDetId(ix, iy - 1, iz)
0261                                           : (EEDetId::validDetId(ix + 1, iy + 1, iz)
0262                                                  ? EEDetId(ix + 1, iy + 1, iz)
0263                                                  : (EEDetId::validDetId(ix + 1, iy - 1, iz)
0264                                                         ? EEDetId(ix + 1, iy - 1, iz)
0265                                                         : (EEDetId::validDetId(ix - 1, iy + 1, iz)
0266                                                                ? EEDetId(ix - 1, iy + 1, iz)
0267                                                                : (EEDetId::validDetId(ix - 1, iy - 1, iz)
0268                                                                       ? EEDetId(ix - 1, iy - 1, iz)
0269                                                                       : did)))))))));
0270 }
0271 
0272 void CaloGeometryAnalyzer::ovrTst(const CaloGeometry* cg,
0273                                   const CaloSubdetectorGeometry* geom,
0274                                   const EEDetId& id,
0275                                   std::fstream& fOvr) {
0276   static const GlobalPoint origin(0, 0, 0);
0277   const int iphi(id.iPhiOuterRing());
0278   if (iphi != 0) {
0279     fOvr << "Barrel Neighbors of Endcap id = " << id << std::endl;
0280     const EcalEndcapGeometry* eeG(dynamic_cast<const EcalEndcapGeometry*>(geom));
0281     auto cell(geom->getGeometry(id));
0282     const CaloSubdetectorGeometry* bar(cg->getSubdetectorGeometry(DetId::Ecal, EcalBarrel));
0283     const EcalEndcapGeometry::OrderedListOfEBDetId* ol(eeG->getClosestBarrelCells(id));
0284     // assert ( nullptr != ol ) ;
0285     if (ol != nullptr) {
0286       for (unsigned int i(0); i != ol->size(); ++i) {
0287         fOvr << "           " << i << "  " << (*ol)[i];
0288         auto other(bar->getGeometry((*ol)[i]));
0289         const GlobalVector cv(cell->getPosition() - origin);
0290         const GlobalVector ov(other->getPosition() - origin);
0291         const double cosang(cv.dot(ov) / (cv.mag() * ov.mag()));
0292         const double angle(convertRadToDeg(acos(std::abs(cosang) < 1. ? cosang : 1.)));
0293         fOvr << ", angle = " << angle << std::endl;
0294       }
0295     }
0296   }
0297 }
0298 
0299 void CaloGeometryAnalyzer::ovrTst(const CaloGeometry* cg,
0300                                   const CaloSubdetectorGeometry* geom,
0301                                   const EBDetId& id,
0302                                   std::fstream& fOvr) {
0303   static const GlobalPoint origin(0, 0, 0);
0304   const int ieta(id.ieta());
0305   if (85 == std::abs(ieta)) {
0306     const EcalBarrelGeometry* ebG(dynamic_cast<const EcalBarrelGeometry*>(geom));
0307     auto cell(geom->getGeometry(id));
0308     const CaloSubdetectorGeometry* ecap(cg->getSubdetectorGeometry(DetId::Ecal, EcalEndcap));
0309     fOvr << "Endcap Neighbors of Barrel id = " << id << std::endl;
0310     const EcalBarrelGeometry::OrderedListOfEEDetId* ol(ebG->getClosestEndcapCells(id));
0311     // assert ( nullptr != ol ) ;
0312     if (ol != nullptr) {
0313       for (unsigned int i(0); i != ol->size(); ++i) {
0314         fOvr << "           " << i << "  " << (*ol)[i];
0315         auto other(ecap->getGeometry((*ol)[i]));
0316         const GlobalVector cv(cell->getPosition() - origin);
0317         const GlobalVector ov(other->getPosition() - origin);
0318         const double cosang(cv.dot(ov) / (cv.mag() * ov.mag()));
0319         const double angle(convertRadToDeg(acos(std::abs(cosang) < 1. ? cosang : 1.)));
0320         fOvr << ", angle = " << angle << std::endl;
0321       }
0322     } else {
0323       fOvr << "endcap ecal ptr is null " << std::endl;
0324     }
0325   }
0326 }
0327 
0328 void CaloGeometryAnalyzer::checkDiff(int i1, int i2, int i3, CenterOrCorner iCtrCor, XorYorZ iXYZ, double diff) {
0329   if (3.5 < fabs(diff)) {
0330     std::cout << "For a volume " << (kCenter == iCtrCor ? "CENTER" : "CORNER") << ", & "
0331               << "i1=" << i1 << " & i2=" << i2 << " & i3=" << i3 << ", ***BIG DISAGREEMENT FOUND. D"
0332               << (kX == iXYZ ? "X" : (kY == iXYZ ? "Y" : "Z")) << "=" << diff << " microns" << std::endl;
0333     m_allOK = false;
0334   }
0335 }
0336 
0337 void CaloGeometryAnalyzer::ctrcor(const DetId& did,
0338                                   const CaloCellGeometry& cell,
0339                                   std::fstream& fCtr,
0340                                   std::fstream& fCor,
0341                                   std::fstream& oldCtr,
0342                                   std::fstream& oldCor,
0343                                   unsigned int histi) {
0344   int oldie(0);
0345   int oldip(0);
0346   oldCtr >> oldie >> oldip;
0347   oldCor >> oldie >> oldip;
0348   const CaloGenericDetId cgid(did);
0349   if (cgid.isEB()) {
0350     const EBDetId ebid(did);
0351     const int ie(ebid.ieta());
0352     const int ip(ebid.iphi());
0353     fCtr << std::setw(4) << ie << std::setw(4) << ip;
0354     fCor << std::setw(4) << ie << std::setw(4) << ip;
0355   }
0356   if (cgid.isEE()) {
0357     const EEDetId eeid(did);
0358     const int ix(eeid.ix());
0359     const int iy(eeid.iy());
0360 
0361     fCtr << std::setw(4) << ix << std::setw(4) << iy;
0362     fCor << std::setw(4) << ix << std::setw(4) << iy;
0363   }
0364   if (cgid.isES()) {
0365     const ESDetId esid(did);
0366     const int pl(esid.plane());
0367     const int ix(esid.six());
0368     const int iy(esid.siy());
0369     const int st(esid.strip());
0370     fCtr << std::setw(4) << pl << std::setw(4) << ix << std::setw(4) << iy << std::setw(4) << st;
0371     fCor << std::setw(4) << pl << std::setw(4) << ix << std::setw(4) << iy << std::setw(4) << st;
0372     int oldiy, oldst;
0373     oldCtr >> oldiy >> oldst;
0374     oldCor >> oldip >> oldst;
0375   }
0376   int depth = 0;
0377   if (cgid.det() == DetId::Hcal) {
0378     const HcalDetId hcid(did);
0379     const int ie(hcid.ieta());
0380     const int ip(hcid.iphi());
0381     const int de(hcid.depth());
0382     fCtr << std::setw(4) << ie << std::setw(4) << ip << std::setw(4) << de;
0383     fCor << std::setw(4) << ie << std::setw(4) << ip << std::setw(4) << de;
0384     int oldde;
0385     oldCtr >> oldde;
0386     oldCor >> oldde;
0387     depth = de;
0388   }
0389   if (cgid.isZDC()) {
0390     const HcalZDCDetId zcid(did);
0391     const int is(zcid.section());
0392     const int ic(zcid.channel());
0393     fCtr << std::setw(4) << is << std::setw(4) << ic;
0394     fCor << std::setw(4) << is << std::setw(4) << ic;
0395   }
0396   if (cgid.isCastor()) {
0397     const HcalCastorDetId cid(did);
0398     const int is(cid.sector());
0399     const int im(cid.module());
0400     fCtr << std::setw(4) << is << std::setw(4) << im;
0401     fCor << std::setw(4) << is << std::setw(4) << im;
0402   }
0403   if (cgid.isCaloTower()) {
0404     const CaloTowerDetId cid(did);
0405     const int ie(cid.ieta());
0406     const int ip(cid.iphi());
0407     fCtr << std::setw(4) << ie << std::setw(4) << ip;
0408     fCor << std::setw(4) << ie << std::setw(4) << ip;
0409   }
0410 
0411   const double x(cell.getPosition().x());
0412   const double y(cell.getPosition().y());
0413   const double z(cell.getPosition().z());
0414 
0415   double oldx, oldy, oldz;
0416 
0417   oldCtr >> oldx >> oldy >> oldz;
0418 
0419   const double dx(1.e4 * (x - oldx));
0420   const double dy(1.e4 * (y - oldy));
0421   const double dz(1.e4 * (z - oldz));
0422 
0423   h_diffs[histi][0]->Fill(dx);
0424   h_diffs[histi][1]->Fill(dy);
0425   h_diffs[histi][2]->Fill(dz);
0426 
0427   checkDiff(oldie, oldip, depth, kCenter, kX, dx);
0428   checkDiff(oldie, oldip, depth, kCenter, kY, dy);
0429   checkDiff(oldie, oldip, depth, kCenter, kZ, dz);
0430 
0431   fCtr << std::fixed << std::setw(12) << std::setprecision(4) << x << std::fixed << std::setw(12)
0432        << std::setprecision(4) << y << std::fixed << std::setw(12) << std::setprecision(4) << z << std::endl;
0433 
0434   const CaloCellGeometry::CornersVec& co(cell.getCorners());
0435 
0436   for (unsigned int j(0); j < co.size(); ++(++(++j))) {
0437     const double x(co[j].x());
0438     const double y(co[j].y());
0439     const double z(co[j].z());
0440 
0441     double oldx, oldy, oldz;
0442 
0443     oldCor >> oldx >> oldy >> oldz;
0444 
0445     const double dx(1.e4 * (x - oldx));
0446     const double dy(1.e4 * (y - oldy));
0447     const double dz(1.e4 * (z - oldz));
0448 
0449     h_diffs[histi][j + 3]->Fill(dx);
0450     h_diffs[histi][j + 4]->Fill(dy);
0451     h_diffs[histi][j + 5]->Fill(dz);
0452 
0453     checkDiff(oldie, oldip, j, kCorner, kX, dx);
0454     checkDiff(oldie, oldip, j, kCorner, kY, dy);
0455     checkDiff(oldie, oldip, j, kCorner, kZ, dz);
0456 
0457     fCor << std::fixed << std::setw(12) << std::setprecision(4) << x << std::fixed << std::setw(12)
0458          << std::setprecision(4) << y << std::fixed << std::setw(12) << std::setprecision(4) << z;
0459   }
0460   fCor << std::endl;
0461 }
0462 
0463 void CaloGeometryAnalyzer::buildHcal(const CaloGeometry* cg,
0464                                      const HcalTopology& ht,
0465                                      DetId::Detector det,
0466                                      int subdetn,
0467                                      const char* name,
0468                                      unsigned int histi) {
0469   std::cout << "Now checking detector " << name << std::endl;
0470   const std::string oldnameCtr("old" + std::string(name) + ".ctr");
0471   const std::string oldnameCor("old" + std::string(name) + ".cor");
0472   const std::string fnameCtr(std::string(name) + ".ctr");
0473   const std::string fnameCor(std::string(name) + ".cor");
0474   const std::string fnameOvr(std::string(name) + ".ovr");
0475   const std::string fnameRoot(std::string(name) + ".C");
0476   std::fstream oldCtr(oldnameCtr.c_str(), std::ios_base::in);
0477   std::fstream oldCor(oldnameCor.c_str(), std::ios_base::in);
0478   std::fstream fCtr(fnameCtr.c_str(), std::ios_base::out);
0479   std::fstream fCor(fnameCor.c_str(), std::ios_base::out);
0480   std::fstream fOvr(fnameOvr.c_str(), std::ios_base::out);
0481   std::fstream f(fnameRoot.c_str(), std::ios_base::out);
0482 
0483   const CaloSubdetectorGeometry* geom(cg->getSubdetectorGeometry(det, subdetn));
0484 
0485   f << "{" << std::endl;
0486   f << "  TGeoManager* geoManager = new TGeoManager(\"ROOT\", \"" << name << "\");" << std::endl;
0487   f << "  TGeoMaterial* dummyMaterial = new TGeoMaterial(\"Vacuum\", 0,0,0); " << std::endl;
0488   f << "  TGeoMedium* dummyMedium =  new TGeoMedium(\"Vacuum\",1,dummyMaterial);" << std::endl;
0489   f << "  TGeoVolume* world=geoManager->MakeBox(\"world\",dummyMedium, 8000.0, 8000.0, 14000.0); " << std::endl;
0490   f << "  geoManager->SetTopVolume(world); " << std::endl;
0491   f << "  TGeoVolume* box; " << std::endl;
0492   int n = 0;
0493   const std::vector<DetId>& ids(geom->getValidDetIds(det, subdetn));
0494 
0495   const std::vector<DetId>& ids2(cg->getValidDetIds(det, subdetn));
0496 
0497   if (ids != ids2) {
0498     std::cout << "Methods differ! One gives size " << ids.size() << " and the other gives size " << ids2.size()
0499               << std::endl;
0500   }
0501   assert(ids == ids2);
0502 
0503   for (const auto& i : ids) {
0504     ++n;
0505     auto cell = (geom->getGeometry(i));
0506 
0507     assert(cg->present(i));
0508 
0509     ctrcor(i, *cell, fCtr, fCor, oldCtr, oldCor, histi);
0510 
0511     const DetId id(i);
0512 
0513     const HcalDetId hcId(i);
0514 
0515     const GlobalPoint pos(cell->getPosition());
0516     const double posmag(pos.mag());
0517 
0518     const double disin(DetId::Ecal == det && EcalPreshower == subdetn ? 0.000001 : 0.001);
0519 
0520     const GlobalPoint pointIn(
0521         pos.x() + disin * pos.x() / posmag, pos.y() + disin * pos.y() / posmag, pos.z() + disin * pos.z() / posmag);
0522     const GlobalPoint pointFr(
0523         pos.x() - 0.1 * pos.x() / posmag, pos.y() - 0.1 * pos.y() / posmag, pos.z() - 0.1 * pos.z() / posmag);
0524 
0525     if (cell->inside(pointFr))
0526       std::cout << "Bad outside: " << pointIn << ", " << pointFr << std::endl;
0527     assert(cell->inside(pointIn));
0528     assert(!cell->inside(pointFr));
0529 
0530     const double deltaPhi(geom->deltaPhi(id));
0531 
0532     const double deltaEta(geom->deltaEta(id));
0533 
0534     const unsigned int detIndex(3);
0535     const GlobalPoint ggp(cell->getPosition());
0536 
0537     h_dPhi[detIndex]->Fill(ht.detId2denseId(hcId), deltaPhi);
0538     h_dPhiR[detIndex]->Fill(ggp.perp(), deltaPhi);
0539 
0540     h_dEta[detIndex]->Fill(ht.detId2denseId(hcId), deltaEta);
0541     h_dEtaR[detIndex]->Fill(ggp.perp(), deltaEta);
0542 
0543     const unsigned int i1(HcalGeometry::alignmentTransformIndexLocal(hcId));
0544 
0545     const DetId d1(HcalGeometry::detIdFromLocalAlignmentIndex(i1));
0546 
0547     const unsigned int i2(HcalGeometry::alignmentTransformIndexLocal(d1));
0548 
0549     assert(i1 == i2);
0550 
0551     f << "  // " << HcalDetId(i) << std::endl;
0552 
0553     const GlobalPoint gp(cell->getPosition());
0554     f << "  // Checking getClosestCell for position " << gp << std::endl;
0555 
0556     HcalDetId closestCell(geom->getClosestCell(gp));
0557 
0558     f << "  // Return position is " << closestCell << std::endl;
0559     if (closestCell != HcalDetId(i)) {
0560       const double rr(reco::deltaR(gp.eta(),
0561                                    gp.phi(),
0562                                    geom->getGeometry(closestCell)->getPosition().eta(),
0563                                    geom->getGeometry(closestCell)->getPosition().phi()));
0564       if (rr > 1.e-5)
0565         std::cout << "For " << HcalDetId(i) << " closest is " << closestCell << " HCAL dR=" << rr << std::endl;
0566     }
0567     // test getCells against base class version every so often
0568     if (0 == ht.detId2denseId(closestCell) % 30) {
0569       cmpset(geom, gp, 2._deg);
0570       cmpset(geom, gp, 5._deg);
0571       cmpset(geom, gp, 7._deg);
0572       cmpset(geom, gp, 25._deg);
0573       cmpset(geom, gp, 45._deg);
0574     }
0575 
0576     if (det == DetId::Hcal && subdetn == HcalForward)
0577       f << "  box=geoManager->MakeBox(\"point\",dummyMedium,1.0,1.0,1.0);" << std::endl;
0578     else
0579       f << "  box=geoManager->MakeBox(\"point\",dummyMedium,3.0,3.0,3.0);" << std::endl;
0580     f << "  world->AddNode(box," << n << ",new TGeoHMatrix(TGeoTranslation(" << cell->getPosition().x() << ","
0581       << cell->getPosition().y() << "," << cell->getPosition().z() << ")));" << std::endl;
0582   }
0583 
0584   f << "  geoManager->CloseGeometry();" << std::endl;
0585   f << "world->Voxelize(\"\"); // now the new geometry is valid for tracking, so you can do \n // even raytracing \n "
0586        "//  if (!canvas) { \n    TCanvas* canvas=new TCanvas(\"EvtDisp\",\"EvtDisp\",500,500); \n //  } \n  "
0587        "canvas->Modified(); \n  canvas->Update();      \n  world->Draw(); \n";
0588   f << "}" << std::endl;
0589   f.close();
0590   fCtr.close();
0591   fCor.close();
0592 }
0593 
0594 void CaloGeometryAnalyzer::build(const CaloGeometry* cg,
0595                                  const HcalTopology& ht,
0596                                  DetId::Detector det,
0597                                  int subdetn,
0598                                  const char* name,
0599                                  unsigned int histi) {
0600   std::cout << "Now checking detector " << name << std::endl;
0601 
0602   const std::string oldnameCtr("old" + std::string(name) + ".ctr");
0603   const std::string oldnameCor("old" + std::string(name) + ".cor");
0604   const std::string fnameCtr(std::string(name) + ".ctr");
0605   const std::string fnameCor(std::string(name) + ".cor");
0606   const std::string fnameOvr(std::string(name) + ".ovr");
0607   const std::string fnameRoot(std::string(name) + ".C");
0608   std::fstream oldCtr(oldnameCtr.c_str(), std::ios_base::in);
0609   std::fstream oldCor(oldnameCor.c_str(), std::ios_base::in);
0610   std::fstream fCtr(fnameCtr.c_str(), std::ios_base::out);
0611   std::fstream fCor(fnameCor.c_str(), std::ios_base::out);
0612   std::fstream fOvr(fnameOvr.c_str(), std::ios_base::out);
0613   std::fstream f(fnameRoot.c_str(), std::ios_base::out);
0614 
0615   const CaloSubdetectorGeometry* geom(cg->getSubdetectorGeometry(det, subdetn));
0616 
0617   f << "{" << std::endl;
0618   f << "  TGeoManager* geoManager = new TGeoManager(\"ROOT\", \"" << name << "\");" << std::endl;
0619   f << "  TGeoMaterial* dummyMaterial = new TGeoMaterial(\"Vacuum\", 0,0,0); " << std::endl;
0620   f << "  TGeoMedium* dummyMedium =  new TGeoMedium(\"Vacuum\",1,dummyMaterial);" << std::endl;
0621   f << "  TGeoVolume* world=geoManager->MakeBox(\"world\",dummyMedium, 8000.0, 8000.0, 14000.0); " << std::endl;
0622   f << "  geoManager->SetTopVolume(world); " << std::endl;
0623   f << "  TGeoVolume* box; " << std::endl;
0624   int n = 0;
0625   const std::vector<DetId>& ids(geom->getValidDetIds(det, subdetn));
0626 
0627   const std::vector<DetId>& ids2(cg->getValidDetIds(det, subdetn));
0628 
0629   if (ids != ids2) {
0630     std::cout << "Methods differ! One gives size " << ids.size() << " and the other gives size " << ids2.size()
0631               << std::endl;
0632   }
0633 
0634   assert(ids == ids2);
0635 
0636   for (const auto& i : ids) {
0637     ++n;
0638     auto cell(geom->getGeometry(i));
0639 
0640     assert(cg->present(i));
0641 
0642     ctrcor(i, *cell, fCtr, fCor, oldCtr, oldCor, histi);
0643 
0644     const DetId id(i);
0645 
0646     const CaloGenericDetId cid(id);
0647 
0648     // assert(cid.validDetId());
0649     // This line fails for CaloTower because this method does not support CaloTower
0650 
0651     // assert(CaloGenericDetId(id.det(), id.subdetId(), cid.denseIndex()) == id);
0652     // This line fails for CaloTower because this three-parameter constructor does not support CaloTower
0653 
0654     const GlobalPoint pos(cell->getPosition());
0655     const double posmag(pos.mag());
0656 
0657     const double disin(DetId::Ecal == det && EcalPreshower == subdetn ? 0.000001 : 0.001);
0658 
0659     const GlobalPoint pointIn(
0660         pos.x() + disin * pos.x() / posmag, pos.y() + disin * pos.y() / posmag, pos.z() + disin * pos.z() / posmag);
0661     const GlobalPoint pointFr(
0662         pos.x() - 0.1 * pos.x() / posmag, pos.y() - 0.1 * pos.y() / posmag, pos.z() - 0.1 * pos.z() / posmag);
0663 
0664     if (cell->inside(pointFr))
0665       std::cout << "Bad outside: " << pointIn << ", " << pointFr << std::endl;
0666     assert(cell->inside(pointIn));
0667     assert(!cell->inside(pointFr));
0668 
0669     const double deltaPhi(geom->deltaPhi(id));
0670 
0671     const double deltaEta(geom->deltaEta(id));
0672 
0673     const unsigned int detIndex(
0674         DetId::Ecal == det && EcalBarrel == subdetn
0675             ? 0
0676             : (DetId::Ecal == det && EcalEndcap == subdetn
0677                    ? 1
0678                    : (DetId::Ecal == det && EcalPreshower == subdetn
0679                           ? 2
0680                           : (DetId::Hcal == det
0681                                  ? 3
0682                                  : (DetId::Calo == det && HcalZDCDetId::SubdetectorId == subdetn
0683                                         ? 4
0684                                         : (DetId::Calo == det && HcalCastorDetId::SubdetectorId == subdetn ? 5 : 6))))));
0685 
0686     const CaloGenericDetId cgid(id);
0687 
0688     if (cgid.isCaloTower() == false) {  // CaloGenericDetId::denseIndex() not supported for CaloTower
0689       const GlobalPoint ggp(cell->getPosition());
0690 
0691       h_dPhi[detIndex]->Fill(cgid.denseIndex(), deltaPhi);
0692       h_dPhiR[detIndex]->Fill(ggp.perp(), deltaPhi);
0693 
0694       h_dEta[detIndex]->Fill(cgid.denseIndex(), deltaEta);
0695       h_dEtaR[detIndex]->Fill(ggp.perp(), deltaEta);
0696     }
0697 
0698     if (det == DetId::Ecal) {
0699       if (subdetn == EcalBarrel) {
0700         f << "  // " << EBDetId(i) << std::endl;
0701 
0702         const GlobalPoint gp(cell->getPosition(0.));
0703 
0704         f << "  // Checking getClosestCell for position " << gp << std::endl;
0705 
0706         const EBDetId ebid(id);
0707         for (unsigned int j(0); j != 4; ++j) {
0708           const CaloCellGeometry::CornersVec& corn(cell->getCorners());
0709           h_eta->Fill(ebid.ieta() * 1. + 0.25 * j, corn[j].eta());
0710           if (ebid.ieta() > 0)
0711             h_phi->Fill(ebid.iphi() * 1. + 0.25 * j, corn[j].phi());
0712         }
0713 
0714         EBDetId closestCell(geom->getClosestCell(gp));
0715 
0716         f << "  // Return position is " << closestCell << std::endl;
0717         assert(closestCell == EBDetId(i));
0718         // test getCells against base class version every so often
0719         if (0 == closestCell.hashedIndex() % 100) {
0720           cmpset(geom, gp, 2._deg);
0721           cmpset(geom, gp, 5._deg);
0722           cmpset(geom, gp, 25._deg);
0723           cmpset(geom, gp, 45._deg);
0724         }
0725 
0726         ovrTst(cg, geom, EBDetId(i), fOvr);
0727 
0728         const unsigned int i1(EcalBarrelGeometry::alignmentTransformIndexLocal(ebid));
0729 
0730         const DetId d1(EcalBarrelGeometry::detIdFromLocalAlignmentIndex(i1));
0731 
0732         const unsigned int i2(EcalBarrelGeometry::alignmentTransformIndexLocal(d1));
0733 
0734         assert(i1 == i2);
0735       }
0736       if (subdetn == EcalEndcap) {
0737         const EEDetId did(i);
0738         const int ix(did.ix());
0739         const int iy(did.iy());
0740         const int iz(did.zside());
0741 
0742         const unsigned int i1(EcalEndcapGeometry::alignmentTransformIndexLocal(did));
0743 
0744         const DetId d1(EcalEndcapGeometry::detIdFromLocalAlignmentIndex(i1));
0745 
0746         const unsigned int i2(EcalEndcapGeometry::alignmentTransformIndexLocal(d1));
0747 
0748         assert(i1 == i2);
0749 
0750         f << "  // Checking getClosestCell for position " << cell->getPosition(0.) << std::endl;
0751 
0752         const GlobalPoint gp(cell->getPosition(0.));
0753 
0754         const EEDetId closestCell(geom->getClosestCell(gp));
0755         f << "  // Return position is " << closestCell << std::endl;
0756 
0757         assert(closestCell == did);
0758         // test getCells against base class version every so often
0759         if (0 == closestCell.hashedIndex() % 10) {
0760           cmpset(geom, gp, 2._deg);
0761           cmpset(geom, gp, 5._deg);
0762           cmpset(geom, gp, 25._deg);
0763           cmpset(geom, gp, 45._deg);
0764         }
0765 
0766         const GlobalVector xx(2.5, 0, 0);
0767         const GlobalVector yy(0, 2.5, 0);
0768         const GlobalVector zz(0, 0, 1);
0769         const GlobalPoint pointIn(cell->getPosition(1.));
0770         const GlobalPoint pointFr(cell->getPosition(-1.));
0771         const GlobalPoint pointBk(cell->getPosition(24.));
0772         const GlobalPoint pointXP(cell->getPosition(1.) + xx);
0773         const GlobalPoint pointXM(cell->getPosition(1.) - xx);
0774         const GlobalPoint pointYP(cell->getPosition(1.) + yy);
0775         const GlobalPoint pointYM(cell->getPosition(1.) - yy);
0776         const GlobalPoint pointPP(cell->getPosition(1.) + xx + yy);
0777         const GlobalPoint pointPM(cell->getPosition(1.) + xx - yy);
0778         const GlobalPoint pointMP(cell->getPosition(1.) - xx + yy);
0779         const GlobalPoint pointMM(cell->getPosition(1.) - xx - yy);
0780         const EEDetId didXP(gid(ix + 1, iy, iz, did));
0781         const EEDetId didXM(gid(ix - 1, iy, iz, did));
0782         const EEDetId didYP(gid(ix, iy + 1, iz, did));
0783         const EEDetId didYM(gid(ix, iy - 1, iz, did));
0784         const EEDetId didPP(gid(ix + 1, iy + 1, iz, did));
0785         const EEDetId didPM(gid(ix + 1, iy - 1, iz, did));
0786         const EEDetId didMP(gid(ix - 1, iy + 1, iz, did));
0787         const EEDetId didMM(gid(ix - 1, iy - 1, iz, did));
0788 
0789         assert(cell->inside(pointIn));
0790         assert(!cell->inside(pointFr));
0791         assert(!cell->inside(pointBk));
0792         assert(!cell->inside(pointXP));
0793         assert(!cell->inside(pointXM));
0794         assert(!cell->inside(pointYP));
0795         assert(!cell->inside(pointYM));
0796         assert(!cell->inside(pointPP));
0797         assert(!cell->inside(pointPM));
0798         assert(!cell->inside(pointMP));
0799         assert(!cell->inside(pointMM));
0800 
0801         const EEDetId ccBk(geom->getClosestCell(pointBk));
0802         const EEDetId ccIn(geom->getClosestCell(pointIn));
0803         const EEDetId ccFr(geom->getClosestCell(pointFr));
0804         const EEDetId ccXP(geom->getClosestCell(pointXP));
0805         const EEDetId ccXM(geom->getClosestCell(pointXM));
0806         const EEDetId ccYP(geom->getClosestCell(pointYP));
0807         const EEDetId ccYM(geom->getClosestCell(pointYM));
0808         const EEDetId ccPP(geom->getClosestCell(pointPP));
0809         const EEDetId ccPM(geom->getClosestCell(pointPM));
0810         const EEDetId ccMP(geom->getClosestCell(pointMP));
0811         const EEDetId ccMM(geom->getClosestCell(pointMM));
0812 
0813         assert(ccIn == did);
0814         assert(ccFr == did);
0815         assert(ccBk == did);
0816         assert(ccXP == didXP || !geom->getGeometry(didXP)->inside(pointXP));
0817         assert(ccXM == didXM || !geom->getGeometry(didXM)->inside(pointXM));
0818         assert(ccYP == didYP || !geom->getGeometry(didYP)->inside(pointYP));
0819         assert(ccYM == didYM || !geom->getGeometry(didYM)->inside(pointYM));
0820         assert(ccPP == didPP || !geom->getGeometry(didPP)->inside(pointPP));
0821         assert(ccPM == didPM || !geom->getGeometry(didPM)->inside(pointPM));
0822         assert(ccMP == didMP || !geom->getGeometry(didMP)->inside(pointMP));
0823         assert(ccMM == didMM || !geom->getGeometry(didMM)->inside(pointMM));
0824 
0825         ovrTst(cg, geom, EEDetId(i), fOvr);
0826       }
0827       if (subdetn == EcalPreshower) {
0828         const ESDetId esid(i);
0829 
0830         f << "  // " << esid << std::endl;
0831         f << "  // Checking getClosestCell for position " << cell->getPosition() << " in plane " << esid.plane()
0832           << std::endl;
0833         ESDetId closestCell = ESDetId((dynamic_cast<const EcalPreshowerGeometry*>(geom))
0834                                           ->getClosestCellInPlane(cell->getPosition(), esid.plane()));
0835         f << "  // Return position is " << closestCell << std::endl;
0836         //sanity checks
0837         int o_zside = esid.zside();
0838 
0839         assert((o_zside < 0 && cell->getPosition().z() < 0.) || (o_zside > 0 && cell->getPosition().z() > 0.));
0840 
0841         if (closestCell != esid)
0842           std::cout << "** esid=" << esid << ", closest=" << closestCell << std::endl;
0843 
0844         const unsigned int i1(EcalPreshowerGeometry::alignmentTransformIndexLocal(esid));
0845 
0846         const DetId d1(EcalPreshowerGeometry::detIdFromLocalAlignmentIndex(i1));
0847 
0848         const unsigned int i2(EcalPreshowerGeometry::alignmentTransformIndexLocal(d1));
0849 
0850         assert(i1 == i2);
0851       }
0852     } else if (det == DetId::Hcal) {
0853       const HcalDetId hcId(i);
0854 
0855       const unsigned int i1(HcalGeometry::alignmentTransformIndexLocal(hcId));
0856 
0857       const DetId d1(HcalGeometry::detIdFromLocalAlignmentIndex(i1));
0858 
0859       const unsigned int i2(HcalGeometry::alignmentTransformIndexLocal(d1));
0860 
0861       assert(i1 == i2);
0862 
0863       f << "  // " << HcalDetId(i) << std::endl;
0864 
0865       const GlobalPoint& gp(cell->getPosition());
0866 
0867       f << "  // Checking getClosestCell for position " << gp << std::endl;
0868 
0869       const HcalDetId closestCell(geom->getClosestCell(gp));
0870 
0871       f << "  // Return position is " << closestCell << std::endl;
0872       if (closestCell != HcalDetId(i)) {
0873         const double rr(reco::deltaR(gp.eta(),
0874                                      gp.phi(),
0875                                      geom->getGeometry(closestCell)->getPosition().eta(),
0876                                      geom->getGeometry(closestCell)->getPosition().phi()));
0877         if (rr > 1.e-5)
0878           std::cout << "For " << HcalDetId(i) << " closest is " << closestCell << " HCAL dR=" << rr << std::endl;
0879       }
0880       // test getCells against base class version every so often
0881       if (0 == ht.detId2denseId(closestCell) % 30) {
0882         cmpset(geom, gp, 2._deg);
0883         cmpset(geom, gp, 5._deg);
0884         cmpset(geom, gp, 7._deg);
0885         cmpset(geom, gp, 25._deg);
0886         cmpset(geom, gp, 45._deg);
0887       }
0888     } else if (det == DetId::Calo && subdetn == HcalCastorDetId::SubdetectorId) {
0889       f << "  // " << HcalCastorDetId(i) << std::endl;
0890 
0891       const GlobalPoint gp(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z() - 0.1);
0892 
0893       f << "  // Checking getClosestCell for position " << gp << std::endl;
0894 
0895       const DetId closestCell(geom->getClosestCell(gp));
0896 
0897       if (closestCell != DetId(0)) {
0898         f << "  // Return position is " << HcalCastorDetId(closestCell) << std::endl;
0899         if (closestCell != HcalCastorDetId(i)) {
0900           const double rr(reco::deltaR(gp.eta(),
0901                                        gp.phi(),
0902                                        geom->getGeometry(closestCell)->getPosition().eta(),
0903                                        geom->getGeometry(closestCell)->getPosition().phi()));
0904           if (rr > 1.e-5)
0905             std::cout << "For " << HcalCastorDetId(i) << " closest is " << HcalCastorDetId(closestCell) << " dR=" << rr
0906                       << std::endl;
0907         }
0908       }
0909       // test getCells against base class version every so often
0910       // if( 0 == closestCell.denseIndex()%30 )
0911       {
0912         cmpset(geom, gp, 2._deg);
0913         cmpset(geom, gp, 5._deg);
0914         cmpset(geom, gp, 7._deg);
0915         cmpset(geom, gp, 25._deg);
0916         cmpset(geom, gp, 45._deg);
0917       }
0918     } else if (det == DetId::Calo && subdetn == HcalZDCDetId::SubdetectorId) {
0919       f << "  // " << HcalZDCDetId(i) << std::endl;
0920       const double sign(HcalZDCDetId(i).zside());
0921       const GlobalPoint gp(cell->getPosition().x(), cell->getPosition().y(), cell->getPosition().z() + sign * 0.1);
0922 
0923       f << "  // Checking getClosestCell for position " << gp << std::endl;
0924 
0925       const DetId closestCell(geom->getClosestCell(gp));
0926 
0927       if (closestCell != DetId(0)) {
0928         f << "  // Return position is " << HcalZDCDetId(closestCell) << std::endl;
0929         if (closestCell != HcalZDCDetId(i)) {
0930           const double rr(reco::deltaR(gp.eta(),
0931                                        gp.phi(),
0932                                        geom->getGeometry(closestCell)->getPosition().eta(),
0933                                        geom->getGeometry(closestCell)->getPosition().phi()));
0934           if (rr > 1.e-5)
0935             std::cout << "For " << HcalZDCDetId(i) << " closest is " << HcalZDCDetId(closestCell) << " dR=" << rr
0936                       << std::endl;
0937         }
0938       }
0939       // test getCells against base class version every so often
0940       // if( 0 == closestCell.denseIndex()%30 )
0941       {
0942         cmpset(geom, gp, 2._deg);
0943         cmpset(geom, gp, 5._deg);
0944         cmpset(geom, gp, 7._deg);
0945         cmpset(geom, gp, 25._deg);
0946         cmpset(geom, gp, 45._deg);
0947       }
0948     }
0949 
0950     if (det == DetId::Hcal && subdetn == HcalForward)
0951       f << "  box=geoManager->MakeBox(\"point\",dummyMedium,1.0,1.0,1.0);" << std::endl;
0952     else
0953       f << "  box=geoManager->MakeBox(\"point\",dummyMedium,3.0,3.0,3.0);" << std::endl;
0954     f << "  world->AddNode(box," << n << ",new TGeoHMatrix(TGeoTranslation(" << cell->getPosition().x() << ","
0955       << cell->getPosition().y() << "," << cell->getPosition().z() << ")));" << std::endl;
0956   }
0957   f << "  geoManager->CloseGeometry();" << std::endl;
0958   f << "world->Voxelize(\"\"); // now the new geometry is valid for tracking, so you can do \n // even raytracing \n "
0959        "//  if (!canvas) { \n    TCanvas* canvas=new TCanvas(\"EvtDisp\",\"EvtDisp\",500,500); \n //  } \n  "
0960        "canvas->Modified(); \n  canvas->Update();      \n  world->Draw(); \n";
0961   f << "}" << std::endl;
0962   f.close();
0963   fCtr.close();
0964   fCor.close();
0965 }
0966 
0967 void CaloGeometryAnalyzer::analyze(const edm::Event& /*iEvent*/, const edm::EventSetup& iSetup) {
0968   const auto& pG = iSetup.getData(geometryToken_);
0969   const CaloGeometry* cG = &pG;
0970   const auto& pT = iSetup.getData(topologyToken_);
0971 
0972   const std::vector<DetId> allDetId(pG.getValidDetIds());
0973 
0974   const std::vector<DetId>& deb(pG.getValidDetIds(DetId::Ecal, EcalBarrel));
0975   const std::vector<DetId>& dee(pG.getValidDetIds(DetId::Ecal, EcalEndcap));
0976   const std::vector<DetId>& des(pG.getValidDetIds(DetId::Ecal, EcalPreshower));
0977   const std::vector<DetId>& dhb(pG.getValidDetIds(DetId::Hcal, HcalBarrel));
0978   const std::vector<DetId>& dhe(pG.getValidDetIds(DetId::Hcal, HcalEndcap));
0979   const std::vector<DetId>& dho(pG.getValidDetIds(DetId::Hcal, HcalOuter));
0980   const std::vector<DetId>& dhf(pG.getValidDetIds(DetId::Hcal, HcalForward));
0981   const std::vector<DetId>& dct(pG.getValidDetIds(DetId::Calo, CaloTowerDetId::SubdetId));
0982   const std::vector<DetId>& dca(pG.getValidDetIds(DetId::Calo, HcalCastorDetId::SubdetectorId));
0983   const std::vector<DetId>& dzd(pG.getValidDetIds(DetId::Calo, HcalZDCDetId::SubdetectorId));
0984 
0985   const std::vector<DetId>& dha(pG.getSubdetectorGeometry(DetId::Hcal, 1)->getValidDetIds());
0986 
0987   const unsigned int sum(deb.size() + dee.size() + des.size() + dhb.size() + dhe.size() + dho.size() + dhf.size() +
0988                          dct.size() + dca.size() + dzd.size());
0989 
0990   if (sum != allDetId.size()) {
0991     std::cout << "Sums differ! One is " << allDetId.size() << " and the other is " << sum << std::endl;
0992   }
0993 
0994   assert(sum == allDetId.size());
0995 
0996   assert(dha.size() == dhb.size() + dhe.size() + dho.size() + dhf.size());
0997   //
0998   // get the ecal & hcal geometry
0999   //
1000   if (pass_ == 0) {
1001     std::cout << "**Ecal Barrel avg Radius = "
1002               << dynamic_cast<const EcalBarrelGeometry*>(pG.getSubdetectorGeometry(DetId::Ecal, EcalBarrel))
1003                      ->avgRadiusXYFrontFaceCenter()
1004               << std::endl;
1005 
1006     std::cout << "**Ecal Endcap avg Zabs = "
1007               << dynamic_cast<const EcalEndcapGeometry*>(pG.getSubdetectorGeometry(DetId::Ecal, EcalEndcap))
1008                      ->avgAbsZFrontFaceCenter()
1009               << std::endl;
1010 
1011     m_allOK = true;
1012 
1013     build(cG, pT, DetId::Ecal, EcalBarrel, "eb", 0);
1014     build(cG, pT, DetId::Ecal, EcalEndcap, "ee", 1);
1015     build(cG, pT, DetId::Ecal, EcalPreshower, "es", 2);
1016     buildHcal(cG, pT, DetId::Hcal, HcalBarrel, "hb", 3);
1017     buildHcal(cG, pT, DetId::Hcal, HcalEndcap, "he", 4);
1018     buildHcal(cG, pT, DetId::Hcal, HcalOuter, "ho", 5);
1019     buildHcal(cG, pT, DetId::Hcal, HcalForward, "hf", 6);
1020     build(cG, pT, DetId::Calo, CaloTowerDetId::SubdetId, "ct", 7);
1021     // build(cG, pT, DetId::Calo, HcalCastorDetId::SubdetectorId, "ca", 8);
1022     // Castor has been removed from CMS
1023     build(cG, pT, DetId::Calo, HcalZDCDetId::SubdetectorId, "zd", 9);
1024 
1025     std::cout << "\n\n*********** Validation of cell centers and corners " << (m_allOK ? "SUCCEEDS!! " : "FAILS!! ")
1026               << "**********************\n\n\n"
1027               << std::endl;
1028   }
1029 
1030   pass_++;
1031 }
1032 
1033 DEFINE_FWK_MODULE(CaloGeometryAnalyzer);