Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2022-05-26 22:39:32 

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