File indexing completed on 2023-03-17 13:02:28
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/MessageLogger/interface/MessageLogger.h"
0029 #include "FWCore/ServiceRegistry/interface/Service.h"
0030 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0031
0032 #include <fstream>
0033 #include <iomanip>
0034 #include <iterator>
0035 #include <sstream>
0036 #include "TH1.h"
0037 #include "TH1D.h"
0038 #include "TProfile.h"
0039
0040 using namespace geant_units;
0041 using namespace geant_units::operators;
0042
0043 class CaloGeometryAnalyzer : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0044 enum CenterOrCorner { kCenter, kCorner };
0045 enum XorYorZ { kX, kY, kZ };
0046
0047 public:
0048 explicit CaloGeometryAnalyzer(const edm::ParameterSet&);
0049 ~CaloGeometryAnalyzer() override;
0050
0051 void beginJob() override {}
0052 void analyze(edm::Event const& iEvent, edm::EventSetup const&) override;
0053 void endJob() override {}
0054
0055 private:
0056
0057 void build(const CaloGeometry* cg,
0058 const HcalTopology& ht,
0059 DetId::Detector det,
0060 int subdetn,
0061 const char* name,
0062 unsigned int histi);
0063
0064 void buildHcal(const CaloGeometry* cg,
0065 const HcalTopology& ht,
0066 DetId::Detector det,
0067 int subdetn,
0068 const char* name,
0069 unsigned int histi);
0070
0071 void ctrcor(const DetId& did,
0072 const CaloCellGeometry& cell,
0073 std::fstream& fCtr,
0074 std::fstream& fCor,
0075 std::fstream& oldCtr,
0076 std::fstream& oldCor,
0077 unsigned int histi);
0078
0079 void checkDiff(int i1, int i2, int i3, CenterOrCorner iCtrCor, XorYorZ iXYZ, double diff);
0080 int pass_;
0081
0082 EEDetId gid(unsigned int ix, unsigned int iy, unsigned int iz, const EEDetId& did) const;
0083
0084 void cmpset(const CaloSubdetectorGeometry* geom, const GlobalPoint& gp, const double dR);
0085
0086 void ovrTst(const CaloGeometry* cg, const CaloSubdetectorGeometry* geom, const EEDetId& id, std::fstream& fOvr);
0087
0088 void ovrTst(const CaloGeometry* cg, const CaloSubdetectorGeometry* geom, const EBDetId& id, std::fstream& fOvr);
0089
0090 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> geometryToken_;
0091 edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> topologyToken_;
0092
0093 edm::Service<TFileService> h_fs;
0094
0095 TProfile* h_dPhi[7];
0096 TProfile* h_dPhiR[7];
0097
0098 TProfile* h_dEta[7];
0099 TProfile* h_dEtaR[7];
0100
0101 TProfile* h_eta;
0102 TProfile* h_phi;
0103
0104 TH1D* h_diffs[10][12];
0105
0106 TH1D* h_scindex;
0107
0108 bool m_allOK;
0109 };
0110
0111 CaloGeometryAnalyzer::CaloGeometryAnalyzer(const edm::ParameterSet& )
0112 : geometryToken_{esConsumes<CaloGeometry, CaloGeometryRecord>(edm::ESInputTag{})},
0113 topologyToken_{esConsumes<HcalTopology, HcalRecNumberingRecord>(edm::ESInputTag{})} {
0114 usesResource("TFileService");
0115
0116 pass_ = 0;
0117
0118 h_dPhi[0] = h_fs->make<TProfile>("dPhi:EB:index", "EB: dPhi vs index", 61200, -0.5, 61199.5, " ");
0119 h_dPhiR[0] = h_fs->make<TProfile>("dPhi:EB:R", "EB: dPhi vs R", 100, 125, 135, " ");
0120
0121 h_dEta[0] = h_fs->make<TProfile>("dEta:EB:index", "EB: dEta vs index", 61200, -0.5, 61199.5, " ");
0122 h_dEtaR[0] = h_fs->make<TProfile>("dEta:EB:R", "EB: dEta vs R", 100, 125, 135, " ");
0123
0124 h_dPhi[1] = h_fs->make<TProfile>("dPhi:EE:index", "EE: dPhi vs index", 14648, -0.5, 14647.5, " ");
0125 h_dPhiR[1] = h_fs->make<TProfile>("dPhi:EE:R", "EE: dPhi vs R", 130, 30, 160, " ");
0126
0127 h_dEta[1] = h_fs->make<TProfile>("dEta:EE:index", "EE: dEta vs index", 14648, -0.5, 14647.5, " ");
0128 h_dEtaR[1] = h_fs->make<TProfile>("dEta:EE:R", "EE: dEta vs R", 130, 30, 160, " ");
0129
0130 h_dPhi[2] = h_fs->make<TProfile>("dPhi:ES:index", "ES: dPhi vs index", 137216, -0.5, 137215.5, " ");
0131 h_dPhiR[2] = h_fs->make<TProfile>("dPhi:ES:R", "ES: dPhi vs R", 90, 40, 130, " ");
0132
0133 h_dEta[2] = h_fs->make<TProfile>("dEta:ES:index", "ES: dEta vs index", 137216, -0.5, 137215.5, " ");
0134 h_dEtaR[2] = h_fs->make<TProfile>("dEta:ES:R", "ES: dEta vs R", 90, 40, 130, " ");
0135
0136 h_dPhi[3] = h_fs->make<TProfile>("dPhi:HC:index", "HC: dPhi vs index", 9072, -0.5, 9071.5, " ");
0137 h_dPhiR[3] = h_fs->make<TProfile>("dPhi:HC:R", "HC: dPhi vs R", 400, 0, 400, " ");
0138
0139 h_dEta[3] = h_fs->make<TProfile>("dEta:HC:index", "HC: dEta vs index", 9072, -0.5, 9071.5, " ");
0140 h_dEtaR[3] = h_fs->make<TProfile>("dEta:HC:R", "HC: dEta vs R", 400, 0, 400, " ");
0141
0142 h_dPhi[4] = h_fs->make<TProfile>("dPhi:ZD:index", "ZD: dPhi vs index", 22, -0.5, 21.5, " ");
0143 h_dPhiR[4] = h_fs->make<TProfile>("dPhi:ZD:R", "ZD: dPhi vs R", 100, 0, 10, " ");
0144
0145 h_dEta[4] = h_fs->make<TProfile>("dEta:ZD:index", "ZD: dEta vs index", 22, -0.5, 21.5, " ");
0146 h_dEtaR[4] = h_fs->make<TProfile>("dEta:ZD:R", "ZD: dEta vs R", 100, 0, 10, " ");
0147
0148 h_dPhi[5] = h_fs->make<TProfile>("dPhi:CA:index", "CA: dPhi vs index", 224, -0.5, 223.5, " ");
0149 h_dPhiR[5] = h_fs->make<TProfile>("dPhi:CA:R", "CA: dPhi vs R", 100, 0, 20, " ");
0150
0151 h_dEta[5] = h_fs->make<TProfile>("dEta:CA:index", "CA: dEta vs index", 224, -0.5, 223.5, " ");
0152 h_dEtaR[5] = h_fs->make<TProfile>("dEta:CA:R", "CA: dEta vs R", 100, 0, 20, " ");
0153
0154 h_dPhi[6] = h_fs->make<TProfile>("dPhi:CT:index", "CT: dPhi vs index", 4320, -0.5, 4319.5, " ");
0155 h_dPhiR[6] = h_fs->make<TProfile>("dPhi:CT:R", "CT: dPhi vs R", 150, 0, 150, " ");
0156
0157 h_dEta[6] = h_fs->make<TProfile>("dEta:CT:index", "CT: dEta vs index", 4320, -0.5, 4319.5, " ");
0158 h_dEtaR[6] = h_fs->make<TProfile>("dEta:CT:R", "CT: dEta vs R", 150, 0, 150, " ");
0159
0160 h_eta = h_fs->make<TProfile>("iEta", "Eta vs iEta", 86 * 2 * 4, -86, 86, " ");
0161 h_phi = h_fs->make<TProfile>("iPhi", "Phi vs iPhi", 360 * 4, 1, 361, " ");
0162
0163 const std::string hname[10] = {"EB", "EE", "ES", "HB", "HO", "HE", "HF", "CT", "ZD", "CA"};
0164 const std::string cname[12] = {
0165 "XCtr", "YCtr", "ZCtr", "XCor0", "YCor0", "ZCor0", "XCor3", "YCor3", "ZCor3", "XCor6", "YCor6", "ZCor6"};
0166
0167 for (unsigned int i(0); i != 10; ++i) {
0168 for (unsigned int j(0); j != 12; ++j) {
0169 h_diffs[i][j] =
0170 h_fs->make<TH1D>(std::string(hname[i] + cname[j] + std::string("Diff (microns)")).c_str(),
0171 std::string(hname[i] + std::string(": New-Nom(") + cname[j] + std::string(")")).c_str(),
0172 200,
0173 -200.,
0174 200.);
0175 }
0176 }
0177 h_scindex = h_fs->make<TH1D>(
0178 std::string("Supercrystal Hashed Index").c_str(), std::string("SC Hashed Index").c_str(), 632, -0.5, 631.5);
0179 }
0180
0181 CaloGeometryAnalyzer::~CaloGeometryAnalyzer() {}
0182
0183 void CaloGeometryAnalyzer::cmpset(const CaloSubdetectorGeometry* geom, const GlobalPoint& gp, const double dR) {
0184 typedef CaloSubdetectorGeometry::DetIdSet DetSet;
0185
0186 DetSet base = geom->CaloSubdetectorGeometry::getCells(gp, dR);
0187 DetSet over = geom->getCells(gp, dR);
0188 if (over == base) {
0189 #ifdef EDM_ML_DEBUG
0190 edm::LogVerbatim("CaloGeom") << "getCells Test dR=" << dR << ", gp=" << gp << ", gp.eta=" << gp.eta()
0191 << ", gp.phi=" << gp.phi()
0192 << ": base and over are equal!\n ***************************\n ";
0193 #endif
0194 } else {
0195 if (2 < std::abs((int)(base.size()) - (int)(over.size()))) {
0196 DetSet inBaseNotOver;
0197 DetSet inOverNotBase;
0198 std::set_difference(
0199 base.begin(), base.end(), over.begin(), over.end(), std::inserter(inBaseNotOver, inBaseNotOver.begin()));
0200
0201 if (inBaseNotOver.empty()) {
0202 edm::LogVerbatim("CaloGeom") << "getCells Test dR=" << dR << ", gp=" << gp << ", gp.eta=" << gp.eta()
0203 << ", gp.phi=" << gp.phi() << ": No elements in base but not overload ";
0204 } else {
0205 edm::LogVerbatim("CaloGeom") << "Length of Base is " << base.size();
0206 edm::LogVerbatim("CaloGeom") << "Length of Over is " << over.size();
0207 edm::LogVerbatim("CaloGeom") << "There are " << inBaseNotOver.size() << " items in Base but not in Overload";
0208
0209 for (const auto& iS : inBaseNotOver) {
0210 std::ostringstream st1;
0211 st1 << "getCells Test dR=" << dR << ", gp=" << gp << ": cell in base but not overload = ";
0212 if (iS.det() == DetId::Ecal && iS.subdetId() == EcalBarrel)
0213 st1 << EBDetId(iS);
0214 if (iS.det() == DetId::Ecal && iS.subdetId() == EcalEndcap)
0215 st1 << EEDetId(iS);
0216 if (iS.det() == DetId::Hcal)
0217 st1 << HcalDetId(iS);
0218 edm::LogVerbatim("CaloGeom") << st1.str();
0219 }
0220 }
0221
0222 std::set_difference(
0223 over.begin(), over.end(), base.begin(), base.end(), std::inserter(inOverNotBase, inOverNotBase.begin()));
0224
0225 if (inOverNotBase.empty()) {
0226 edm::LogVerbatim("CaloGeom") << "getCells Test dR=" << dR << ", gp=" << gp << ", gp.eta=" << gp.eta()
0227 << ", gp.phi=" << gp.phi() << ": No elements in overload but not base ";
0228 } else {
0229 edm::LogVerbatim("CaloGeom") << "Length of Base is " << base.size();
0230 edm::LogVerbatim("CaloGeom") << "Length of Over is " << over.size();
0231 edm::LogVerbatim("CaloGeom") << "There are " << inOverNotBase.size() << " items in Overload but not in Base";
0232
0233 for (const auto& iS : inOverNotBase) {
0234 std::ostringstream st1;
0235 st1 << "getCells Test dR=" << dR << ", gp=" << gp << ": cell in overload but not base = ";
0236 if (iS.det() == DetId::Ecal && iS.subdetId() == EcalBarrel)
0237 st1 << EBDetId(iS);
0238 if (iS.det() == DetId::Ecal && iS.subdetId() == EcalEndcap)
0239 st1 << EEDetId(iS);
0240 if (iS.det() == DetId::Hcal)
0241 st1 << HcalDetId(iS);
0242 edm::LogVerbatim("CaloGeom") << st1.str();
0243 }
0244 }
0245 edm::LogVerbatim("CaloGeom") << "------------- done with mismatch printout ---------------";
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
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
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 edm::LogVerbatim("CaloGeom") << "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";
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 edm::LogVerbatim("CaloGeom") << "Now checking detector " << name;
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 edm::LogVerbatim("CaloGeom") << "Methods differ! One gives size " << ids.size() << " and the other gives size "
0499 << ids2.size();
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 edm::LogVerbatim("CaloGeom") << "Bad outside: " << pointIn << ", " << pointFr;
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 edm::LogVerbatim("CaloGeom") << "For " << HcalDetId(i) << " closest is " << closestCell << " HCAL dR=" << rr;
0566 }
0567
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 edm::LogVerbatim("CaloGeom") << "Now checking detector " << name;
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 edm::LogVerbatim("CaloGeom") << "Methods differ! One gives size " << ids.size() << " and the other gives size "
0631 << ids2.size();
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
0649
0650
0651
0652
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 edm::LogVerbatim("CaloGeom") << "Bad outside: " << pointIn << ", " << pointFr;
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) {
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
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
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
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 edm::LogVerbatim("CaloGeom") << "** esid=" << esid << ", closest=" << closestCell;
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 edm::LogVerbatim("CaloGeom") << "For " << HcalDetId(i) << " closest is " << closestCell << " HCAL dR=" << rr;
0879 }
0880
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 edm::LogVerbatim("CaloGeom") << "For " << HcalCastorDetId(i) << " closest is "
0906 << HcalCastorDetId(closestCell) << " dR=" << rr;
0907 }
0908 }
0909
0910
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 edm::LogVerbatim("CaloGeom") << "For " << HcalZDCDetId(i) << " closest is " << HcalZDCDetId(closestCell)
0936 << " dR=" << rr;
0937 }
0938 }
0939
0940
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& , 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 edm::LogVerbatim("CaloGeom") << "Sums differ! One is " << allDetId.size() << " and the other is " << sum;
0992 }
0993
0994 assert(sum == allDetId.size());
0995
0996 assert(dha.size() == dhb.size() + dhe.size() + dho.size() + dhf.size());
0997
0998
0999
1000 if (pass_ == 0) {
1001 edm::LogVerbatim("CaloGeom") << "**Ecal Barrel avg Radius = "
1002 << dynamic_cast<const EcalBarrelGeometry*>(
1003 pG.getSubdetectorGeometry(DetId::Ecal, EcalBarrel))
1004 ->avgRadiusXYFrontFaceCenter();
1005
1006 edm::LogVerbatim("CaloGeom") << "**Ecal Endcap avg Zabs = "
1007 << dynamic_cast<const EcalEndcapGeometry*>(
1008 pG.getSubdetectorGeometry(DetId::Ecal, EcalEndcap))
1009 ->avgAbsZFrontFaceCenter();
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
1022
1023 build(cG, pT, DetId::Calo, HcalZDCDetId::SubdetectorId, "zd", 9);
1024
1025 edm::LogVerbatim("CaloGeom") << "\n\n*********** Validation of cell centers and corners "
1026 << (m_allOK ? "SUCCEEDS!! " : "FAILS!! ") << "**********************\n\n\n";
1027 }
1028
1029 pass_++;
1030 }
1031
1032 DEFINE_FWK_MODULE(CaloGeometryAnalyzer);
