Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-11-02 23:48:53 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!