Geometry/src/FWRecoGeometryESProducer.cc

0001 #include <fstream>
0002 #include <streambuf>
0003
0004 #include "Fireworks/Geometry/interface/FWRecoGeometryESProducer.h"
0005 #include "Fireworks/Geometry/interface/FWRecoGeometry.h"
0006 #include "Fireworks/Geometry/interface/FWTGeoRecoGeometry.h"
0007 #include "Fireworks/Geometry/interface/FWRecoGeometryRecord.h"
0008
0009 #include "DataFormats/GeometrySurface/interface/RectangularPlaneBounds.h"
0010 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0011 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0012 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0013 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0014 #include "Geometry/HGCalGeometry/interface/HGCalGeometry.h"
0015 #include "Geometry/MTDGeometryBuilder/interface/MTDGeometry.h"
0016 #include "Geometry/MTDGeometryBuilder/interface/ProxyMTDTopology.h"
0017 #include "Geometry/MTDGeometryBuilder/interface/RectangularMTDTopology.h"
0018 #include "RecoLocalCalo/HGCalRecAlgos/interface/RecHitTools.h"
0019 #include "Geometry/CaloGeometry/interface/CaloCellGeometry.h"
0020 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0021 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0022 #include "Geometry/CSCGeometry/interface/CSCChamber.h"
0023 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0024 #include "Geometry/DTGeometry/interface/DTChamber.h"
0025 #include "Geometry/DTGeometry/interface/DTLayer.h"
0026 #include "Geometry/RPCGeometry/interface/RPCGeometry.h"
0027 #include "Geometry/GEMGeometry/interface/GEMGeometry.h"
0028 #include "Geometry/GEMGeometry/interface/ME0Geometry.h"
0029 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0030 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0031 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0032 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0033 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0034 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetType.h"
0035 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0036 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0037 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
0038 #include "Geometry/CommonTopologies/interface/TrapezoidalStripTopology.h"
0039 #include "Geometry/CommonTopologies/interface/GEMStripTopology.h"
0040
0041 #include "TNamed.h"
0042 #include "FWCore/ParameterSet/interface/FileInPath.h"
0043
0044 void FWRecoGeometryESProducer::ADD_PIXEL_TOPOLOGY(unsigned int rawid,
0045                                                   const GeomDet* detUnit,
0046                                                   FWRecoGeometry& fwRecoGeometry) {
0047   const PixelGeomDetUnit* det = dynamic_cast<const PixelGeomDetUnit*>(detUnit);
0048   if (det) {
0049     const PixelTopology* topo = &det->specificTopology();
0050
0051     std::pair<float, float> pitch = topo->pitch();
0052     fwRecoGeometry.idToName[rawid].topology[0] = pitch.first;
0053     fwRecoGeometry.idToName[rawid].topology[1] = pitch.second;
0054
0055     fwRecoGeometry.idToName[rawid].topology[2] = topo->localX(0.f);  // offsetX
0056     fwRecoGeometry.idToName[rawid].topology[3] = topo->localY(0.f);  // offsetY
0057
0058     // big pixels layout
0059     fwRecoGeometry.idToName[rawid].topology[4] = topo->isItBigPixelInX(80) ? 0 : 1;
0060   }
0061 }
0062
0063 using Phase2TrackerGeomDetUnit = PixelGeomDetUnit;
0064 using Phase2TrackerTopology = PixelTopology;
0065
0066 #define ADD_SISTRIP_TOPOLOGY(rawid, detUnit)                                                                   \
0067   const StripGeomDetUnit* det = dynamic_cast<const StripGeomDetUnit*>(detUnit);                                \
0068   if (det) {                                                                                                   \
0069     if (const StripTopology* topo = dynamic_cast<const StripTopology*>(&det->specificTopology())) {            \
0070       fwRecoGeometry.idToName[rawid].topology[0] = 0;                                                          \
0071       fwRecoGeometry.idToName[rawid].topology[1] = topo->nstrips();                                            \
0072       fwRecoGeometry.idToName[rawid].topology[2] = topo->stripLength();                                        \
0073     }                                                                                                          \
0074     if (const RadialStripTopology* rtop =                                                                      \
0075             dynamic_cast<const RadialStripTopology*>(&(det->specificType().specificTopology()))) {             \
0076       fwRecoGeometry.idToName[rawid].topology[0] = 1;                                                          \
0077       fwRecoGeometry.idToName[rawid].topology[3] = rtop->yAxisOrientation();                                   \
0078       fwRecoGeometry.idToName[rawid].topology[4] = rtop->originToIntersection();                               \
0079       fwRecoGeometry.idToName[rawid].topology[5] = rtop->phiOfOneEdge();                                       \
0080       fwRecoGeometry.idToName[rawid].topology[6] = rtop->angularWidth();                                       \
0081     } else if (const RectangularStripTopology* topo =                                                          \
0082                    dynamic_cast<const RectangularStripTopology*>(&(det->specificType().specificTopology()))) { \
0083       fwRecoGeometry.idToName[rawid].topology[0] = 2;                                                          \
0084       fwRecoGeometry.idToName[rawid].topology[3] = topo->pitch();                                              \
0085     } else if (const TrapezoidalStripTopology* topo =                                                          \
0086                    dynamic_cast<const TrapezoidalStripTopology*>(&(det->specificType().specificTopology()))) { \
0087       fwRecoGeometry.idToName[rawid].topology[0] = 3;                                                          \
0088       fwRecoGeometry.idToName[rawid].topology[3] = topo->pitch();                                              \
0089     }                                                                                                          \
0090   } else {                                                                                                     \
0091     const Phase2TrackerGeomDetUnit* det = dynamic_cast<const Phase2TrackerGeomDetUnit*>(detUnit);              \
0092     if (det) {                                                                                                 \
0093       if (const Phase2TrackerTopology* topo =                                                                  \
0094               dynamic_cast<const Phase2TrackerTopology*>(&(det->specificTopology()))) {                        \
0095         fwRecoGeometry.idToName[rawid].topology[0] = topo->pitch().first;                                      \
0096         fwRecoGeometry.idToName[rawid].topology[1] = topo->pitch().second;                                     \
0097       }                                                                                                        \
0098     }                                                                                                          \
0099   }
0100
0101 void FWRecoGeometryESProducer::ADD_MTD_TOPOLOGY(unsigned int rawid,
0102                                                 const GeomDet* detUnit,
0103                                                 FWRecoGeometry& fwRecoGeometry) {
0104   const MTDGeomDet* det = dynamic_cast<const MTDGeomDet*>(detUnit);
0105
0106   if (det) {
0107     const ProxyMTDTopology& topoproxy = static_cast<const ProxyMTDTopology&>(det->topology());
0108     const RectangularMTDTopology& topo = static_cast<const RectangularMTDTopology&>(topoproxy.specificTopology());
0109
0110     std::pair<float, float> pitch = topo.pitch();
0111     fwRecoGeometry.idToName[rawid].topology[0] = pitch.first;
0112     fwRecoGeometry.idToName[rawid].topology[1] = pitch.second;
0113
0114     fwRecoGeometry.idToName[rawid].topology[2] = topo.xoffset();
0115     fwRecoGeometry.idToName[rawid].topology[3] = topo.yoffset();
0116   }
0117 }
0118
0119 namespace {
0120   const std::array<std::string, 3> hgcal_geom_names = {
0121       {"HGCalEESensitive", "HGCalHESiliconSensitive", "HGCalHEScintillatorSensitive"}};
0122 }
0123
0124 FWRecoGeometryESProducer::FWRecoGeometryESProducer(const edm::ParameterSet& pset) : m_current(-1) {
0125   m_tracker = pset.getUntrackedParameter<bool>("Tracker", true);
0126   m_muon = pset.getUntrackedParameter<bool>("Muon", true);
0127   m_gem = pset.getUntrackedParameter<bool>("GEM", false);
0128   m_calo = pset.getUntrackedParameter<bool>("Calo", true);
0129   m_timing = pset.getUntrackedParameter<bool>("Timing", false);
0130   auto cc = setWhatProduced(this);
0131
0132   if (m_muon)
0133     m_gem = true;
0134   if (m_tracker or m_muon or m_gem) {
0135     m_trackingGeomToken = cc.consumes();
0136   }
0137   if (m_calo) {
0138     m_caloGeomToken = cc.consumes();
0139   }
0140   if (m_timing) {
0141     m_mtdGeomToken = cc.consumes();
0142   }
0143 }
0144
0145 FWRecoGeometryESProducer::~FWRecoGeometryESProducer(void) {}
0146
0147 std::unique_ptr<FWRecoGeometry> FWRecoGeometryESProducer::produce(const FWRecoGeometryRecord& record) {
0148   using namespace edm;
0149
0150   auto fwRecoGeometry = std::make_unique<FWRecoGeometry>();
0151
0152   if (m_tracker || m_muon || m_gem) {
0153     m_trackingGeom = &record.get(m_trackingGeomToken);
0154   }
0155
0156   if (m_tracker) {
0157     DetId detId(DetId::Tracker, 0);
0158     m_trackerGeom = static_cast<const TrackerGeometry*>(m_trackingGeom->slaveGeometry(detId));
0159     addPixelBarrelGeometry(*fwRecoGeometry);
0160     addPixelForwardGeometry(*fwRecoGeometry);
0161     addTIBGeometry(*fwRecoGeometry);
0162     addTIDGeometry(*fwRecoGeometry);
0163     addTOBGeometry(*fwRecoGeometry);
0164     addTECGeometry(*fwRecoGeometry);
0165     writeTrackerParametersXML(*fwRecoGeometry);
0166   }
0167   if (m_muon) {
0168     addDTGeometry(*fwRecoGeometry);
0169     addCSCGeometry(*fwRecoGeometry);
0170     addRPCGeometry(*fwRecoGeometry);
0171     addME0Geometry(*fwRecoGeometry);
0172   }
0173   if (m_gem) {
0174     addGEMGeometry(*fwRecoGeometry);
0175   }
0176   if (m_calo) {
0177     m_caloGeom = &record.get(m_caloGeomToken);
0178     addCaloGeometry(*fwRecoGeometry);
0179   }
0180   if (m_timing) {
0181     m_mtdGeom = &record.get(m_mtdGeomToken);
0182     addMTDGeometry(*fwRecoGeometry);
0183   }
0184
0185   fwRecoGeometry->idToName.resize(m_current + 1);
0186   std::vector<FWRecoGeom::Info>(fwRecoGeometry->idToName).swap(fwRecoGeometry->idToName);
0187   std::sort(fwRecoGeometry->idToName.begin(), fwRecoGeometry->idToName.end());
0188
0189   return fwRecoGeometry;
0190 }
0191
0192 void FWRecoGeometryESProducer::addCSCGeometry(FWRecoGeometry& fwRecoGeometry) {
0193   DetId detId(DetId::Muon, 2);
0194   const CSCGeometry* cscGeometry = static_cast<const CSCGeometry*>(m_trackingGeom->slaveGeometry(detId));
0195   for (auto it = cscGeometry->chambers().begin(), end = cscGeometry->chambers().end(); it != end; ++it) {
0196     const CSCChamber* chamber = *it;
0197
0198     if (chamber) {
0199       unsigned int rawid = chamber->geographicalId();
0200       unsigned int current = insert_id(rawid, fwRecoGeometry);
0201       fillShapeAndPlacement(current, chamber, fwRecoGeometry);
0202       //
0203       // CSC layers geometry
0204       //
0205       for (std::vector<const CSCLayer*>::const_iterator lit = chamber->layers().begin(), lend = chamber->layers().end();
0206            lit != lend;
0207            ++lit) {
0208         const CSCLayer* layer = *lit;
0209
0210         if (layer) {
0211           unsigned int rawid = layer->geographicalId();
0212           unsigned int current = insert_id(rawid, fwRecoGeometry);
0213           fillShapeAndPlacement(current, layer, fwRecoGeometry);
0214
0215           const CSCStripTopology* stripTopology = layer->geometry()->topology();
0216           fwRecoGeometry.idToName[current].topology[0] = stripTopology->yAxisOrientation();
0217           fwRecoGeometry.idToName[current].topology[1] = stripTopology->centreToIntersection();
0218           fwRecoGeometry.idToName[current].topology[2] = stripTopology->yCentreOfStripPlane();
0219           fwRecoGeometry.idToName[current].topology[3] = stripTopology->phiOfOneEdge();
0220           fwRecoGeometry.idToName[current].topology[4] = stripTopology->stripOffset();
0221           fwRecoGeometry.idToName[current].topology[5] = stripTopology->angularWidth();
0222
0223           const CSCWireTopology* wireTopology = layer->geometry()->wireTopology();
0224           fwRecoGeometry.idToName[current].topology[6] = wireTopology->wireSpacing();
0225           fwRecoGeometry.idToName[current].topology[7] = wireTopology->wireAngle();
0226         }
0227       }
0228     }
0229   }
0230 }
0231
0232 void FWRecoGeometryESProducer::addDTGeometry(FWRecoGeometry& fwRecoGeometry) {
0233   DetId detId(DetId::Muon, 1);
0234   const DTGeometry* dtGeometry = static_cast<const DTGeometry*>(m_trackingGeom->slaveGeometry(detId));
0235
0236   //
0237   // DT chambers geometry
0238   //
0239   for (auto it = dtGeometry->chambers().begin(), end = dtGeometry->chambers().end(); it != end; ++it) {
0240     const DTChamber* chamber = *it;
0241
0242     if (chamber) {
0243       unsigned int rawid = chamber->geographicalId().rawId();
0244       unsigned int current = insert_id(rawid, fwRecoGeometry);
0245       fillShapeAndPlacement(current, chamber, fwRecoGeometry);
0246     }
0247   }
0248
0249   // Fill in DT layer parameters
0250   for (auto it = dtGeometry->layers().begin(), end = dtGeometry->layers().end(); it != end; ++it) {
0251     const DTLayer* layer = *it;
0252
0253     if (layer) {
0254       unsigned int rawid = layer->id().rawId();
0255       unsigned int current = insert_id(rawid, fwRecoGeometry);
0256       fillShapeAndPlacement(current, layer, fwRecoGeometry);
0257
0258       const DTTopology& topo = layer->specificTopology();
0259       const BoundPlane& surf = layer->surface();
0260       // Topology W/H/L:
0261       fwRecoGeometry.idToName[current].topology[0] = topo.cellWidth();
0262       fwRecoGeometry.idToName[current].topology[1] = topo.cellHeight();
0263       fwRecoGeometry.idToName[current].topology[2] = topo.cellLenght();
0264       fwRecoGeometry.idToName[current].topology[3] = topo.firstChannel();
0265       fwRecoGeometry.idToName[current].topology[4] = topo.lastChannel();
0266       fwRecoGeometry.idToName[current].topology[5] = topo.channels();
0267
0268       // Bounds W/H/L:
0269       fwRecoGeometry.idToName[current].topology[6] = surf.bounds().width();
0270       fwRecoGeometry.idToName[current].topology[7] = surf.bounds().thickness();
0271       fwRecoGeometry.idToName[current].topology[8] = surf.bounds().length();
0272     }
0273   }
0274 }
0275
0276 void FWRecoGeometryESProducer::addRPCGeometry(FWRecoGeometry& fwRecoGeometry) {
0277   //
0278   // RPC rolls geometry
0279   //
0280   DetId detId(DetId::Muon, 3);
0281   const RPCGeometry* rpcGeom = static_cast<const RPCGeometry*>(m_trackingGeom->slaveGeometry(detId));
0282   for (auto it = rpcGeom->rolls().begin(), end = rpcGeom->rolls().end(); it != end; ++it) {
0283     const RPCRoll* roll = (*it);
0284     if (roll) {
0285       unsigned int rawid = roll->geographicalId().rawId();
0286       unsigned int current = insert_id(rawid, fwRecoGeometry);
0287       fillShapeAndPlacement(current, roll, fwRecoGeometry);
0288
0289       const StripTopology& topo = roll->specificTopology();
0290       fwRecoGeometry.idToName[current].topology[0] = topo.nstrips();
0291       fwRecoGeometry.idToName[current].topology[1] = topo.stripLength();
0292       fwRecoGeometry.idToName[current].topology[2] = topo.pitch();
0293     }
0294   }
0295
0296   try {
0297     RPCDetId id(1, 1, 4, 1, 1, 1, 1);
0298     m_trackingGeom->slaveGeometry(detId);
0299     fwRecoGeometry.extraDet.Add(new TNamed("RE4", "RPC endcap station 4"));
0300   } catch (std::runtime_error& e) {
0301     std::cerr << e.what() << std::endl;
0302   }
0303 }
0304
0305 void FWRecoGeometryESProducer::addGEMGeometry(FWRecoGeometry& fwRecoGeometry) {
0306   //
0307   // GEM geometry
0308   //
0309
0310   try {
0311     DetId detId(DetId::Muon, 4);
0312     const GEMGeometry* gemGeom = static_cast<const GEMGeometry*>(m_trackingGeom->slaveGeometry(detId));
0313
0314     // add in superChambers - gem Segments are based on superChambers
0315     for (auto sc : gemGeom->superChambers()) {
0316       if (sc) {
0317         unsigned int rawid = sc->geographicalId().rawId();
0318         unsigned int current = insert_id(rawid, fwRecoGeometry);
0319         fillShapeAndPlacement(current, sc, fwRecoGeometry);
0320       }
0321     }
0322     // add in chambers
0323     for (auto ch : gemGeom->chambers()) {
0324       if (ch) {
0325         unsigned int rawid = ch->geographicalId().rawId();
0326         unsigned int current = insert_id(rawid, fwRecoGeometry);
0327         fillShapeAndPlacement(current, ch, fwRecoGeometry);
0328       }
0329     }
0330     // add in etaPartitions - gem rechits are based on etaPartitions
0331     for (auto roll : gemGeom->etaPartitions()) {
0332       if (roll) {
0333         unsigned int rawid = roll->geographicalId().rawId();
0334         unsigned int current = insert_id(rawid, fwRecoGeometry);
0335         fillShapeAndPlacement(current, roll, fwRecoGeometry);
0336
0337         const StripTopology& topo = roll->specificTopology();
0338         fwRecoGeometry.idToName[current].topology[0] = topo.nstrips();
0339         fwRecoGeometry.idToName[current].topology[1] = topo.stripLength();
0340         fwRecoGeometry.idToName[current].topology[2] = topo.pitch();
0341
0342         float height = topo.stripLength() / 2;
0343         LocalPoint lTop(0., height, 0.);
0344         LocalPoint lBottom(0., -height, 0.);
0345         fwRecoGeometry.idToName[current].topology[3] = roll->localPitch(lTop);
0346         fwRecoGeometry.idToName[current].topology[4] = roll->localPitch(lBottom);
0347         fwRecoGeometry.idToName[current].topology[5] = roll->npads();
0348       }
0349     }
0350
0351     fwRecoGeometry.extraDet.Add(new TNamed("GEM", "GEM muon detector"));
0352     try {
0353       GEMDetId id(1, 1, 2, 1, 1, 1);
0354       m_trackingGeom->slaveGeometry(id);
0355       fwRecoGeometry.extraDet.Add(new TNamed("GE2", "GEM endcap station 2"));
0356     } catch (std::runtime_error& e) {
0357       std::cerr << e.what() << std::endl;
0358     }
0359
0360   } catch (cms::Exception& exception) {
0361     edm::LogError("FWRecoGeometry") << " GEM geometry not found " << exception.what() << std::endl;
0362   }
0363 }
0364
0365 void FWRecoGeometryESProducer::addME0Geometry(FWRecoGeometry& fwRecoGeometry) {
0366   //
0367   // ME0 geometry
0368   //
0369
0370   DetId detId(DetId::Muon, 5);
0371   try {
0372     const ME0Geometry* me0Geom = static_cast<const ME0Geometry*>(m_trackingGeom->slaveGeometry(detId));
0373     for (auto roll : me0Geom->etaPartitions()) {
0374       if (roll) {
0375         unsigned int rawid = roll->geographicalId().rawId();
0376         unsigned int current = insert_id(rawid, fwRecoGeometry);
0377         fillShapeAndPlacement(current, roll, fwRecoGeometry);
0378
0379         const StripTopology& topo = roll->specificTopology();
0380         fwRecoGeometry.idToName[current].topology[0] = topo.nstrips();
0381         fwRecoGeometry.idToName[current].topology[1] = topo.stripLength();
0382         fwRecoGeometry.idToName[current].topology[2] = topo.pitch();
0383
0384         float height = topo.stripLength() / 2;
0385         LocalPoint lTop(0., height, 0.);
0386         LocalPoint lBottom(0., -height, 0.);
0387         fwRecoGeometry.idToName[current].topology[3] = roll->localPitch(lTop);
0388         fwRecoGeometry.idToName[current].topology[4] = roll->localPitch(lBottom);
0389         fwRecoGeometry.idToName[current].topology[5] = roll->npads();
0390       }
0391     }
0392     fwRecoGeometry.extraDet.Add(new TNamed("ME0", "ME0 muon detector"));
0393   } catch (cms::Exception& exception) {
0394     edm::LogError("FWRecoGeometry") << " ME0 geometry not found " << exception.what() << std::endl;
0395   }
0396 }
0397
0398 void FWRecoGeometryESProducer::addPixelBarrelGeometry(FWRecoGeometry& fwRecoGeometry) {
0399   for (TrackerGeometry::DetContainer::const_iterator it = m_trackerGeom->detsPXB().begin(),
0400                                                      end = m_trackerGeom->detsPXB().end();
0401        it != end;
0402        ++it) {
0403     const GeomDet* det = *it;
0404
0405     if (det) {
0406       DetId detid = det->geographicalId();
0407       unsigned int rawid = detid.rawId();
0408       unsigned int current = insert_id(rawid, fwRecoGeometry);
0409       fillShapeAndPlacement(current, det, fwRecoGeometry);
0410
0411       ADD_PIXEL_TOPOLOGY(current, m_trackerGeom->idToDetUnit(detid), fwRecoGeometry);
0412     }
0413   }
0414 }
0415
0416 void FWRecoGeometryESProducer::addPixelForwardGeometry(FWRecoGeometry& fwRecoGeometry) {
0417   for (TrackerGeometry::DetContainer::const_iterator it = m_trackerGeom->detsPXF().begin(),
0418                                                      end = m_trackerGeom->detsPXF().end();
0419        it != end;
0420        ++it) {
0421     const GeomDet* det = *it;
0422
0423     if (det) {
0424       DetId detid = det->geographicalId();
0425       unsigned int rawid = detid.rawId();
0426       unsigned int current = insert_id(rawid, fwRecoGeometry);
0427       fillShapeAndPlacement(current, det, fwRecoGeometry);
0428
0429       ADD_PIXEL_TOPOLOGY(current, m_trackerGeom->idToDetUnit(detid), fwRecoGeometry);
0430     }
0431   }
0432 }
0433
0434 void FWRecoGeometryESProducer::addTIBGeometry(FWRecoGeometry& fwRecoGeometry) {
0435   for (TrackerGeometry::DetContainer::const_iterator it = m_trackerGeom->detsTIB().begin(),
0436                                                      end = m_trackerGeom->detsTIB().end();
0437        it != end;
0438        ++it) {
0439     const GeomDet* det = *it;
0440
0441     if (det) {
0442       DetId detid = det->geographicalId();
0443       unsigned int rawid = detid.rawId();
0444       unsigned int current = insert_id(rawid, fwRecoGeometry);
0445       fillShapeAndPlacement(current, det, fwRecoGeometry);
0446
0447       ADD_SISTRIP_TOPOLOGY(current, m_trackerGeom->idToDet(detid));
0448     }
0449   }
0450 }
0451
0452 void FWRecoGeometryESProducer::addTOBGeometry(FWRecoGeometry& fwRecoGeometry) {
0453   for (TrackerGeometry::DetContainer::const_iterator it = m_trackerGeom->detsTOB().begin(),
0454                                                      end = m_trackerGeom->detsTOB().end();
0455        it != end;
0456        ++it) {
0457     const GeomDet* det = *it;
0458
0459     if (det) {
0460       DetId detid = det->geographicalId();
0461       unsigned int rawid = detid.rawId();
0462       unsigned int current = insert_id(rawid, fwRecoGeometry);
0463       fillShapeAndPlacement(current, det, fwRecoGeometry);
0464
0465       ADD_SISTRIP_TOPOLOGY(current, m_trackerGeom->idToDet(detid));
0466     }
0467   }
0468 }
0469
0470 void FWRecoGeometryESProducer::addTIDGeometry(FWRecoGeometry& fwRecoGeometry) {
0471   for (TrackerGeometry::DetContainer::const_iterator it = m_trackerGeom->detsTID().begin(),
0472                                                      end = m_trackerGeom->detsTID().end();
0473        it != end;
0474        ++it) {
0475     const GeomDet* det = *it;
0476
0477     if (det) {
0478       DetId detid = det->geographicalId();
0479       unsigned int rawid = detid.rawId();
0480       unsigned int current = insert_id(rawid, fwRecoGeometry);
0481       fillShapeAndPlacement(current, det, fwRecoGeometry);
0482
0483       ADD_SISTRIP_TOPOLOGY(current, m_trackerGeom->idToDet(detid));
0484     }
0485   }
0486 }
0487
0488 void FWRecoGeometryESProducer::addTECGeometry(FWRecoGeometry& fwRecoGeometry) {
0489   for (TrackerGeometry::DetContainer::const_iterator it = m_trackerGeom->detsTEC().begin(),
0490                                                      end = m_trackerGeom->detsTEC().end();
0491        it != end;
0492        ++it) {
0493     const GeomDet* det = *it;
0494
0495     if (det) {
0496       DetId detid = det->geographicalId();
0497       unsigned int rawid = detid.rawId();
0498       unsigned int current = insert_id(rawid, fwRecoGeometry);
0499       fillShapeAndPlacement(current, det, fwRecoGeometry);
0500
0501       ADD_SISTRIP_TOPOLOGY(current, m_trackerGeom->idToDet(detid));
0502     }
0503   }
0504 }
0505
0506 void FWRecoGeometryESProducer::addCaloGeometry(FWRecoGeometry& fwRecoGeometry) {
0507   std::vector<DetId> vid = m_caloGeom->getValidDetIds();  // Calo
0508   std::set<DetId> cache;
0509   for (std::vector<DetId>::const_iterator it = vid.begin(), end = vid.end(); it != end; ++it) {
0510     unsigned int id = insert_id(it->rawId(), fwRecoGeometry);
0511     if (!((DetId::Forward == it->det()) || (DetId::HGCalEE == it->det()) || (DetId::HGCalHSi == it->det()) ||
0512           (DetId::HGCalHSc == it->det()))) {
0513       const CaloCellGeometry::CornersVec& cor = m_caloGeom->getGeometry(*it)->getCorners();
0514       fillPoints(id, cor.begin(), cor.end(), fwRecoGeometry);
0515     } else {
0516       DetId::Detector det = it->det();
0517       int subdet = (((DetId::HGCalEE == det) || (DetId::HGCalHSi == det) || (DetId::HGCalHSc == det)) ? ForwardEmpty
0518                                                                                                       : it->subdetId());
0519       const HGCalGeometry* geom = dynamic_cast<const HGCalGeometry*>(m_caloGeom->getSubdetectorGeometry(det, subdet));
0520       hgcal::RecHitTools rhtools;
0521       rhtools.setGeometry(*m_caloGeom);
0522       const auto cor = geom->getNewCorners(*it);
0523
0524       // roll = yaw = pitch = 0
0525       fwRecoGeometry.idToName[id].matrix[0] = 1.0;
0526       fwRecoGeometry.idToName[id].matrix[4] = 1.0;
0527       fwRecoGeometry.idToName[id].matrix[8] = 1.0;
0528
0529       // corners of the front face
0530       for (uint i = 0; i < (cor.size() - 1); ++i) {
0531         fwRecoGeometry.idToName[id].points[i * 3 + 0] = cor[i].x();
0532         fwRecoGeometry.idToName[id].points[i * 3 + 1] = cor[i].y();
0533         fwRecoGeometry.idToName[id].points[i * 3 + 2] = cor[i].z();
0534       }
0535
0536       // center
0537       auto center = geom->getPosition(*it);
0538       fwRecoGeometry.idToName[id].points[(cor.size() - 1) * 3 + 0] = center.x();
0539       fwRecoGeometry.idToName[id].points[(cor.size() - 1) * 3 + 1] = center.y();
0540       fwRecoGeometry.idToName[id].points[(cor.size() - 1) * 3 + 2] = center.z();
0541
0542       // Cells rotation (read few lines below)
0543       fwRecoGeometry.idToName[id].shape[2] = 0.;
0544
0545       // Thickness
0546       fwRecoGeometry.idToName[id].shape[3] = cor[cor.size() - 1].z();
0547
0548       // total points
0549       fwRecoGeometry.idToName[id].topology[0] = cor.size() - 1;
0550
0551       // Layer with Offset
0552       fwRecoGeometry.idToName[id].topology[1] = rhtools.getLayerWithOffset(it->rawId());
0553
0554       // Zside, +/- 1
0555       fwRecoGeometry.idToName[id].topology[2] = rhtools.zside(it->rawId());
0556
0557       // Is Silicon
0558       fwRecoGeometry.idToName[id].topology[3] = rhtools.isSilicon(it->rawId());
0559
0560       // Silicon index
0561       fwRecoGeometry.idToName[id].topology[4] =
0562           rhtools.isSilicon(it->rawId()) ? rhtools.getSiThickIndex(it->rawId()) : -1.;
0563
0564       // Last EE layer
0565       fwRecoGeometry.idToName[id].topology[5] = rhtools.lastLayerEE();
0566
0567       // Compute the orientation of each cell. The orientation here is simply
0568       // addressing the corner or side bottom layout of the cell and should not
0569       // be confused with the concept of orientation embedded in the flat-file
0570       // description. The default orientation of the cells within a wafer is
0571       // with the side at the bottom. The points returned by the HGCal query
0572       // will be ordered counter-clockwise, with the first corner in the
0573       // uppermost-right position. The corners used to calculate the angle wrt
0574       // the Y scale are corner 0 and corner 3, that are opposite in the cells.
0575       // The angle should be 30 degrees wrt the Y axis for all cells in the
0576       // default position. For the rotated layers in CE-H, the situation is
0577       // such that the cells are oriented with a vertex down (assuming those
0578       // layers will have a 30 degrees rotation): this will establish an angle
0579       // of 60 degrees wrt the Y axis. The default way in which an hexagon is
0580       // rendered inside Fireworks is with the vertex down.
0581       if (rhtools.isSilicon(it->rawId())) {
0582         auto val_x = (cor[0].x() - cor[3].x());
0583         auto val_y = (cor[0].y() - cor[3].y());
0584         auto val = round(std::acos(val_y / std::sqrt(val_x * val_x + val_y * val_y)) / M_PI * 180.);
0585         // Pass down the chain the vaue of the rotation of the cell wrt the Y axis.
0586         fwRecoGeometry.idToName[id].shape[2] = val;
0587       }
0588
0589       // For each and every wafer in HGCal, add a "fake" DetId with cells'
0590       // (u,v) bits set to 1 . Those DetIds will be used inside Fireworks to
0591       // render the HGCal Geometry. Due to the huge number of cells involved,
0592       // the HGCal geometry for the Silicon Sensor is wafer-based, not cells
0593       // based. The representation of the single RecHits and of all quantities
0594       // derived from those, is instead fully cells based. The geometry
0595       // representation of the Scintillator is directly based on tiles,
0596       // therefore no fake DetId creations is needed.
0597       if ((det == DetId::HGCalEE) || (det == DetId::HGCalHSi)) {
0598         // Avoid hard coding masks by using static data members from HGCSiliconDetId
0599         auto maskZeroUV = (HGCSiliconDetId::kHGCalCellVMask << HGCSiliconDetId::kHGCalCellVOffset) |
0600                           (HGCSiliconDetId::kHGCalCellUMask << HGCSiliconDetId::kHGCalCellUOffset);
0601         DetId wafer_detid = it->rawId() | maskZeroUV;
0602         // Just be damn sure that's a fake id.
0603         assert(wafer_detid != it->rawId());
0604         auto [iter, is_new] = cache.insert(wafer_detid);
0605         if (is_new) {
0606           unsigned int local_id = insert_id(wafer_detid, fwRecoGeometry);
0607           auto const& dddConstants = geom->topology().dddConstants();
0608           auto wafer_size = static_cast<float>(dddConstants.waferSize(true));
0609           auto R = wafer_size / std::sqrt(3.f);
0610           auto r = wafer_size / 2.f;
0611           float x[6] = {-r, -r, 0.f, r, r, 0.f};
0612           float y[6] = {R / 2.f, -R / 2.f, -R, -R / 2.f, R / 2.f, R};
0613           float z[6] = {0.f, 0.f, 0.f, 0.f, 0.f, 0.f};
0614           for (unsigned int i = 0; i < 6; ++i) {
0615             HepGeom::Point3D<float> wafer_corner(x[i], y[i], z[i]);
0616             auto point =
0617                 geom->topology().dddConstants().waferLocal2Global(wafer_corner, wafer_detid, true, true, false);
0618             fwRecoGeometry.idToName[local_id].points[i * 3 + 0] = point.x();
0619             fwRecoGeometry.idToName[local_id].points[i * 3 + 1] = point.y();
0620             fwRecoGeometry.idToName[local_id].points[i * 3 + 2] = point.z();
0621           }
0622           // Nota Bene: rotations of full layers (and wafers therein) is taken
0623           // care of internally by the call to the waferLocal2Global method.
0624           // Therefore we set up the unit matrix for the rotation.
0625           // roll = yaw = pitch = 0
0626           fwRecoGeometry.idToName[local_id].matrix[0] = 1.0;
0627           fwRecoGeometry.idToName[local_id].matrix[4] = 1.0;
0628           fwRecoGeometry.idToName[local_id].matrix[8] = 1.0;
0629
0630           // thickness
0631           fwRecoGeometry.idToName[local_id].shape[3] = cor[cor.size() - 1].z();
0632
0633           // total points
0634           fwRecoGeometry.idToName[local_id].topology[0] = 6;
0635
0636           // Layer with Offset
0637           fwRecoGeometry.idToName[local_id].topology[1] = rhtools.getLayerWithOffset(it->rawId());
0638
0639           // Zside, +/- 1
0640           fwRecoGeometry.idToName[local_id].topology[2] = rhtools.zside(it->rawId());
0641
0642           // Is Silicon
0643           fwRecoGeometry.idToName[local_id].topology[3] = rhtools.isSilicon(it->rawId());
0644
0645           // Silicon index
0646           fwRecoGeometry.idToName[local_id].topology[4] =
0647               rhtools.isSilicon(it->rawId()) ? rhtools.getSiThickIndex(it->rawId()) : -1.;
0648
0649           // Last EE layer
0650           fwRecoGeometry.idToName[local_id].topology[5] = rhtools.lastLayerEE();
0651         }
0652       }
0653     }
0654   }
0655 }
0656
0657 void FWRecoGeometryESProducer::addMTDGeometry(FWRecoGeometry& fwRecoGeometry) {
0658   for (auto const& det : m_mtdGeom->detUnits()) {
0659     if (det) {
0660       DetId detid = det->geographicalId();
0661       unsigned int rawid = detid.rawId();
0662       unsigned int current = insert_id(rawid, fwRecoGeometry);
0663       fillShapeAndPlacement(current, det, fwRecoGeometry);
0664
0665       ADD_MTD_TOPOLOGY(current, m_mtdGeom->idToDetUnit(detid), fwRecoGeometry);
0666     }
0667   }
0668 }
0669
0670 unsigned int FWRecoGeometryESProducer::insert_id(unsigned int rawid, FWRecoGeometry& fwRecoGeometry) {
0671   ++m_current;
0672   fwRecoGeometry.idToName.push_back(FWRecoGeom::Info());
0673   fwRecoGeometry.idToName.back().id = rawid;
0674
0675   return m_current;
0676 }
0677
0678 void FWRecoGeometryESProducer::fillPoints(unsigned int id,
0679                                           std::vector<GlobalPoint>::const_iterator begin,
0680                                           std::vector<GlobalPoint>::const_iterator end,
0681                                           FWRecoGeometry& fwRecoGeometry) {
0682   unsigned int index(0);
0683   for (std::vector<GlobalPoint>::const_iterator i = begin; i != end; ++i) {
0684     assert(index < FWTGeoRecoGeometry::maxPoints_ - 1);
0685     fwRecoGeometry.idToName[id].points[index] = i->x();
0686     fwRecoGeometry.idToName[id].points[++index] = i->y();
0687     fwRecoGeometry.idToName[id].points[++index] = i->z();
0688     ++index;
0689   }
0690 }
0691
0692 /** Shape of GeomDet */
0693 void FWRecoGeometryESProducer::fillShapeAndPlacement(unsigned int id,
0694                                                      const GeomDet* det,
0695                                                      FWRecoGeometry& fwRecoGeometry) {
0696   // Trapezoidal
0697   const Bounds* b = &((det->surface()).bounds());
0698   if (const TrapezoidalPlaneBounds* b2 = dynamic_cast<const TrapezoidalPlaneBounds*>(b)) {
0699     std::array<const float, 4> const& par = b2->parameters();
0700
0701     // These parameters are half-lengths, as in CMSIM/GEANT3
0702     fwRecoGeometry.idToName[id].shape[0] = 1;
0703     fwRecoGeometry.idToName[id].shape[1] = par[0];  // hBottomEdge
0704     fwRecoGeometry.idToName[id].shape[2] = par[1];  // hTopEdge
0705     fwRecoGeometry.idToName[id].shape[3] = par[2];  // thickness
0706     fwRecoGeometry.idToName[id].shape[4] = par[3];  // apothem
0707   }
0708   if (const RectangularPlaneBounds* b2 = dynamic_cast<const RectangularPlaneBounds*>(b)) {
0709     // Rectangular
0710     fwRecoGeometry.idToName[id].shape[0] = 2;
0711     fwRecoGeometry.idToName[id].shape[1] = b2->width() * 0.5;      // half width
0712     fwRecoGeometry.idToName[id].shape[2] = b2->length() * 0.5;     // half length
0713     fwRecoGeometry.idToName[id].shape[3] = b2->thickness() * 0.5;  // half thickness
0714   }
0715
0716   // Position of the DetUnit's center
0717   GlobalPoint pos = det->surface().position();
0718   fwRecoGeometry.idToName[id].translation[0] = pos.x();
0719   fwRecoGeometry.idToName[id].translation[1] = pos.y();
0720   fwRecoGeometry.idToName[id].translation[2] = pos.z();
0721
0722   // Add the coeff of the rotation matrix
0723   // with a projection on the basis vectors
0724   TkRotation<float> detRot = det->surface().rotation();
0725   fwRecoGeometry.idToName[id].matrix[0] = detRot.xx();
0726   fwRecoGeometry.idToName[id].matrix[1] = detRot.yx();
0727   fwRecoGeometry.idToName[id].matrix[2] = detRot.zx();
0728   fwRecoGeometry.idToName[id].matrix[3] = detRot.xy();
0729   fwRecoGeometry.idToName[id].matrix[4] = detRot.yy();
0730   fwRecoGeometry.idToName[id].matrix[5] = detRot.zy();
0731   fwRecoGeometry.idToName[id].matrix[6] = detRot.xz();
0732   fwRecoGeometry.idToName[id].matrix[7] = detRot.yz();
0733   fwRecoGeometry.idToName[id].matrix[8] = detRot.zz();
0734 }
0735
0736 void FWRecoGeometryESProducer::writeTrackerParametersXML(FWRecoGeometry& fwRecoGeometry) {
0737   std::string path = "Geometry/TrackerCommonData/data/";
0738   if (m_trackerGeom->isThere(GeomDetEnumerators::P1PXB) || m_trackerGeom->isThere(GeomDetEnumerators::P1PXEC)) {
0739     path += "PhaseI/";
0740   } else if (m_trackerGeom->isThere(GeomDetEnumerators::P2PXB) || m_trackerGeom->isThere(GeomDetEnumerators::P2PXEC) ||
0741              m_trackerGeom->isThere(GeomDetEnumerators::P2OTB) || m_trackerGeom->isThere(GeomDetEnumerators::P2OTEC)) {
0742     path += "PhaseII/";
0743   }
0744   path += "trackerParameters.xml";
0745   std::string fullPath = edm::FileInPath(path).fullPath();
0746   std::ifstream t(fullPath);
0747   std::stringstream buffer;
0748   buffer << t.rdbuf();
0749   fwRecoGeometry.trackerTopologyXML = buffer.str();
0750 }