Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-07-12 02:58:35

0001 //#define EDM_ML_DEBUG
0002 
0003 #include "RecoMTD/DetLayers/interface/ETLDetLayerGeometryBuilder.h"
0004 
0005 #include <RecoMTD/DetLayers/interface/MTDRingForwardDoubleLayer.h>
0006 #include <RecoMTD/DetLayers/interface/MTDDetRing.h>
0007 #include <RecoMTD/DetLayers/interface/MTDSectorForwardDoubleLayer.h>
0008 #include <RecoMTD/DetLayers/interface/MTDDetSector.h>
0009 #include <DataFormats/ForwardDetId/interface/ETLDetId.h>
0010 #include <Geometry/CommonDetUnit/interface/GeomDet.h>
0011 #include <Geometry/MTDCommonData/interface/MTDTopologyMode.h>
0012 
0013 #include <Utilities/General/interface/precomputed_value_sort.h>
0014 #include <Geometry/CommonDetUnit/interface/DetSorting.h>
0015 #include <FWCore/MessageLogger/interface/MessageLogger.h>
0016 
0017 #include <iostream>
0018 
0019 using namespace std;
0020 
0021 pair<vector<DetLayer*>, vector<DetLayer*> > ETLDetLayerGeometryBuilder::buildLayers(const MTDGeometry& geo,
0022                                                                                     const MTDTopology& topo) {
0023   vector<DetLayer*> result[2];  // one for each endcap
0024 
0025   const int mtdTopologyMode = topo.getMTDTopologyMode();
0026   ETLDetId::EtlLayout etlL = MTDTopologyMode::etlLayoutFromTopoMode(mtdTopologyMode);
0027   if (etlL == ETLDetId::EtlLayout::tp) {
0028     for (unsigned endcap = 0; endcap < 2; ++endcap) {
0029       // there is only one layer for ETL right now, maybe more later
0030       for (unsigned layer = 0; layer < ETLDetId::kETLv1nDisc; ++layer) {
0031         vector<unsigned> rings;
0032         rings.reserve(ETLDetId::kETLv1maxRing + 1);
0033         for (unsigned ring = 1; ring <= ETLDetId::kETLv1maxRing; ++ring) {
0034           rings.push_back(ring);
0035         }
0036         MTDRingForwardDoubleLayer* thelayer = buildLayer(endcap, layer, rings, geo);
0037         if (thelayer)
0038           result[endcap].push_back(thelayer);
0039       }
0040     }
0041   } else {
0042     // number of layers is identical for post TDR scenarios, pick v4
0043     // loop on number of sectors per face, two faces per disc (i.e. layer) taken into account in layer building (front/back)
0044     unsigned int nSector(1);
0045     if (etlL == ETLDetId::EtlLayout::v4) {
0046       nSector *= ETLDetId::kETLv4maxSector;
0047     } else if (etlL == ETLDetId::EtlLayout::v5) {
0048       nSector *= ETLDetId::kETLv5maxSector;
0049     } else {
0050       throw cms::Exception("MTDDetLayers") << "Not implemented scenario " << mtdTopologyMode;
0051     }
0052 
0053     for (unsigned endcap = 0; endcap < 2; ++endcap) {
0054       // number of layers is two, identical for post TDR scenarios, pick v4
0055       for (unsigned layer = 1; layer <= ETLDetId::kETLv4nDisc; ++layer) {
0056         vector<unsigned> sectors;
0057         sectors.reserve(nSector + 1);
0058         for (unsigned sector = 1; sector <= nSector; ++sector) {
0059           sectors.push_back(sector);
0060         }
0061         MTDSectorForwardDoubleLayer* thelayer = buildLayerNew(endcap, layer, sectors, geo, topo);
0062         if (thelayer)
0063           result[endcap].push_back(thelayer);
0064       }
0065     }
0066   }
0067   //
0068   // the first entry is Z+ ( MTD side 1), the second is Z- (MTD side 0)
0069   //
0070   pair<vector<DetLayer*>, vector<DetLayer*> > res_pair(result[1], result[0]);
0071   return res_pair;
0072 }
0073 
0074 MTDRingForwardDoubleLayer* ETLDetLayerGeometryBuilder::buildLayer(int endcap,
0075                                                                   int layer,
0076                                                                   vector<unsigned>& rings,
0077                                                                   const MTDGeometry& geo) {
0078   MTDRingForwardDoubleLayer* result = nullptr;
0079 
0080   vector<const ForwardDetRing*> frontRings, backRings;
0081 
0082   for (unsigned ring : rings) {
0083     vector<const GeomDet*> frontGeomDets, backGeomDets;
0084     for (unsigned module = 1; module <= ETLDetId::kETLmoduleMask; ++module) {
0085       ETLDetId detId(endcap, ring, module, 0);
0086       const GeomDet* geomDet = geo.idToDet(detId);
0087       // we sometimes loop over more chambers than there are in ring
0088       bool isInFront = isFront(layer, ring, module);
0089       if (geomDet != nullptr) {
0090         if (isInFront) {
0091           frontGeomDets.push_back(geomDet);
0092         } else {
0093           backGeomDets.push_back(geomDet);
0094         }
0095         LogTrace("MTDDetLayers") << "get ETL module " << std::hex << ETLDetId(endcap, layer, ring, module).rawId()
0096                                  << std::dec << " at R=" << geomDet->position().perp()
0097                                  << ", phi=" << geomDet->position().phi() << ", z= " << geomDet->position().z()
0098                                  << " isFront? " << isInFront;
0099       }
0100     }
0101 
0102     if (!backGeomDets.empty()) {
0103       backRings.push_back(makeDetRing(backGeomDets));
0104     }
0105 
0106     if (!frontGeomDets.empty()) {
0107       frontRings.push_back(makeDetRing(frontGeomDets));
0108       assert(!backGeomDets.empty());
0109       float frontz = frontRings[0]->position().z();
0110       float backz = backRings[0]->position().z();
0111       assert(std::abs(frontz) < std::abs(backz));
0112     }
0113   }
0114 
0115   // How should they be sorted?
0116   //    precomputed_value_sort(muDetRods.begin(), muDetRods.end(), geomsort::ExtractZ<GeometricSearchDet,float>());
0117   result = new MTDRingForwardDoubleLayer(frontRings, backRings);
0118   LogTrace("MTDDetLayers") << "New MTDRingForwardLayer with " << frontRings.size() << " and " << backRings.size()
0119                            << " rings, at Z " << result->position().z()
0120                            << " R1: " << result->specificSurface().innerRadius()
0121                            << " R2: " << result->specificSurface().outerRadius();
0122   return result;
0123 }
0124 
0125 bool ETLDetLayerGeometryBuilder::isFront(int layer, int ring, int module) { return (module + 1) % 2; }
0126 
0127 MTDDetRing* ETLDetLayerGeometryBuilder::makeDetRing(vector<const GeomDet*>& geomDets) {
0128   precomputed_value_sort(geomDets.begin(), geomDets.end(), geomsort::DetPhi());
0129   MTDDetRing* result = new MTDDetRing(geomDets);
0130   LogTrace("MTDDetLayers") << "ETLDetLayerGeometryBuilder: new MTDDetRing with " << geomDets.size()
0131                            << " chambers at z=" << result->position().z()
0132                            << " R1: " << result->specificSurface().innerRadius()
0133                            << " R2: " << result->specificSurface().outerRadius();
0134   return result;
0135 }
0136 
0137 MTDSectorForwardDoubleLayer* ETLDetLayerGeometryBuilder::buildLayerNew(
0138     int endcap, int layer, vector<unsigned>& sectors, const MTDGeometry& geo, const MTDTopology& topo) {
0139   MTDSectorForwardDoubleLayer* result = nullptr;
0140 
0141   std::vector<const MTDDetSector*> frontSectors, backSectors;
0142 
0143   LogDebug("MTDDetLayers") << "ETL dets array size = " << geo.detsETL().size();
0144 
0145   for (unsigned sector : sectors) {
0146     std::vector<const GeomDet*> frontGeomDets, backGeomDets;
0147     LogDebug("MTDDetLayers") << "endcap = " << endcap << " layer = " << layer << " sector = " << sector;
0148 #ifdef EDM_ML_DEBUG
0149     unsigned int nfront(0), nback(0);
0150 #endif
0151     for (auto det : geo.detsETL()) {
0152       ETLDetId theMod(det->geographicalId().rawId());
0153       if (theMod.mtdSide() == endcap && theMod.nDisc() == layer && theMod.sector() == static_cast<int>(sector)) {
0154         LogTrace("MTDLayerDump") << std::fixed << "ETLDetId " << theMod.rawId() << " side = " << std::setw(4)
0155                                  << theMod.mtdSide() << " Disc/Side/Sector = " << std::setw(4) << theMod.nDisc() << " "
0156                                  << std::setw(4) << theMod.discSide() << " " << std::setw(4) << theMod.sector()
0157                                  << " mod/type = " << std::setw(4) << theMod.module() << " " << std::setw(4)
0158                                  << theMod.modType() << " pos = " << det->position();
0159         // front layer face
0160         if (theMod.discSide() == 0) {
0161 #ifdef EDM_ML_DEBUG
0162           nfront++;
0163           LogTrace("MTDDetLayers") << "Front " << theMod.discSide() << " " << nfront;
0164 #endif
0165           frontGeomDets.emplace_back(det);
0166           // back layer face
0167         } else if (theMod.discSide() == 1) {
0168 #ifdef EDM_ML_DEBUG
0169           nback++;
0170           LogTrace("MTDDetLayers") << "Back " << theMod.discSide() << " " << nback;
0171 #endif
0172           backGeomDets.emplace_back(det);
0173         }
0174       }
0175     }
0176 
0177     if (!backGeomDets.empty()) {
0178       std::sort(backGeomDets.begin(), backGeomDets.end(), topo.orderETLSector);
0179       LogDebug("MTDDetLayers") << "backGeomDets size = " << backGeomDets.size();
0180       backSectors.emplace_back(makeDetSector(backGeomDets, topo));
0181     }
0182 
0183     if (!frontGeomDets.empty()) {
0184       std::sort(frontGeomDets.begin(), frontGeomDets.end(), topo.orderETLSector);
0185       LogDebug("MTDDetLayers") << "frontGeomDets size = " << frontGeomDets.size();
0186       frontSectors.emplace_back(makeDetSector(frontGeomDets, topo));
0187       assert(!backGeomDets.empty());
0188       float frontz = frontSectors.back()->position().z();
0189       float backz = backSectors.back()->position().z();
0190       assert(std::abs(frontz) < std::abs(backz));
0191     }
0192   }
0193 
0194   result = new MTDSectorForwardDoubleLayer(frontSectors, backSectors);
0195   LogTrace("MTDDetLayers") << "New MTDSectorForwardDoubleLayer with " << std::fixed << std::setw(14)
0196                            << frontSectors.size() << " and " << std::setw(14) << backSectors.size() << " rings, at Z "
0197                            << std::setw(14) << result->specificSurface().position().z() << " R1: " << std::setw(14)
0198                            << result->specificSurface().innerRadius() << " R2: " << std::setw(14)
0199                            << result->specificSurface().outerRadius();
0200 
0201   return result;
0202 }
0203 
0204 MTDDetSector* ETLDetLayerGeometryBuilder::makeDetSector(vector<const GeomDet*>& geomDets, const MTDTopology& topo) {
0205   MTDDetSector* result = new MTDDetSector(geomDets, topo);
0206   LogTrace("MTDDetLayers") << "ETLDetLayerGeometryBuilder::makeDetSector new MTDDetSector with " << std::fixed
0207                            << std::setw(14) << geomDets.size() << " modules \n"
0208                            << (*result);
0209 
0210   return result;
0211 }