LSTCore/src/ModuleMethods.h

0001 #ifndef RecoTracker_LSTCore_src_ModuleMethods_h
0002 #define RecoTracker_LSTCore_src_ModuleMethods_h
0003
0004 #include <map>
0005 #include <iostream>
0006
0007 #include "RecoTracker/LSTCore/interface/Common.h"
0008 #include "RecoTracker/LSTCore/interface/ModulesSoA.h"
0009 #include "RecoTracker/LSTCore/interface/ModulesHostCollection.h"
0010 #include "RecoTracker/LSTCore/interface/TiltedGeometry.h"
0011 #include "RecoTracker/LSTCore/interface/EndcapGeometry.h"
0012 #include "RecoTracker/LSTCore/interface/ModuleConnectionMap.h"
0013 #include "RecoTracker/LSTCore/interface/PixelMap.h"
0014
0015 #include "HeterogeneousCore/AlpakaInterface/interface/host.h"
0016 #include "HeterogeneousCore/AlpakaInterface/interface/memory.h"
0017
0018 namespace lst {
0019   struct ModuleMetaData {
0020     std::map<unsigned int, uint16_t> detIdToIndex;
0021     std::map<unsigned int, float> module_x;
0022     std::map<unsigned int, float> module_y;
0023     std::map<unsigned int, float> module_z;
0024     std::map<unsigned int, unsigned int> module_type;  // 23 : Ph2PSP, 24 : Ph2PSS, 25 : Ph2SS
0025     // https://github.com/cms-sw/cmssw/blob/5e809e8e0a625578aa265dc4b128a93830cb5429/Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h#L29
0026   };
0027
0028   bool parseIsLower(bool isInvertedx, unsigned int detId) { return (isInvertedx) ? !(detId & 1) : (detId & 1); }
0029
0030   unsigned int parsePartnerModuleId(unsigned int detId, bool isLowerx, bool isInvertedx) {
0031     return isLowerx ? (isInvertedx ? detId - 1 : detId + 1) : (isInvertedx ? detId + 1 : detId - 1);
0032   }
0033
0034   bool parseIsInverted(short subdet, short side, short module, short layer) {
0035     if (subdet == Endcap) {
0036       if (side == NegZ) {
0037         return module % 2 == 1;
0038       } else if (side == PosZ) {
0039         return module % 2 == 0;
0040       } else {
0041         return false;
0042       }
0043     } else if (subdet == Barrel) {
0044       if (side == Center) {
0045         if (layer <= 3) {
0046           return module % 2 == 1;
0047         } else if (layer >= 4) {
0048           return module % 2 == 0;
0049         } else {
0050           return false;
0051         }
0052       } else if (side == NegZ or side == PosZ) {
0053         if (layer <= 2) {
0054           return module % 2 == 1;
0055         } else if (layer == 3) {
0056           return module % 2 == 0;
0057         } else {
0058           return false;
0059         }
0060       } else {
0061         return false;
0062       }
0063     } else {
0064       return false;
0065     }
0066   }
0067
0068   inline std::tuple<unsigned int,
0069                     std::vector<unsigned int>,
0070                     unsigned int,
0071                     std::vector<unsigned int>,
0072                     unsigned int,
0073                     std::vector<unsigned int>>
0074   getConnectedPixels(uint16_t nModules, unsigned int& nPixels, PixelMap& pixelMapping, MapPLStoLayer const& pLStoLayer) {
0075     std::vector<unsigned int> connectedModuleDetIds;
0076     std::vector<unsigned int> connectedModuleDetIds_pos;
0077     std::vector<unsigned int> connectedModuleDetIds_neg;
0078
0079     unsigned int totalSizes = 0;
0080     unsigned int totalSizes_pos = 0;
0081     unsigned int totalSizes_neg = 0;
0082     for (unsigned int isuperbin = 0; isuperbin < size_superbins; isuperbin++) {
0083       int sizes = 0;
0084       for (auto const& mCM_pLS : pLStoLayer[0]) {
0085         std::vector<unsigned int> connectedModuleDetIds_pLS =
0086             mCM_pLS.getConnectedModuleDetIds(isuperbin + size_superbins);
0087         connectedModuleDetIds.insert(
0088             connectedModuleDetIds.end(), connectedModuleDetIds_pLS.begin(), connectedModuleDetIds_pLS.end());
0089         sizes += connectedModuleDetIds_pLS.size();
0090       }
0091       pixelMapping.connectedPixelsIndex[isuperbin] = totalSizes;
0092       pixelMapping.connectedPixelsSizes[isuperbin] = sizes;
0093       totalSizes += sizes;
0094
0095       int sizes_pos = 0;
0096       for (auto const& mCM_pLS : pLStoLayer[1]) {
0097         std::vector<unsigned int> connectedModuleDetIds_pLS_pos = mCM_pLS.getConnectedModuleDetIds(isuperbin);
0098         connectedModuleDetIds_pos.insert(connectedModuleDetIds_pos.end(),
0099                                          connectedModuleDetIds_pLS_pos.begin(),
0100                                          connectedModuleDetIds_pLS_pos.end());
0101         sizes_pos += connectedModuleDetIds_pLS_pos.size();
0102       }
0103       pixelMapping.connectedPixelsIndexPos[isuperbin] = totalSizes_pos;
0104       pixelMapping.connectedPixelsSizesPos[isuperbin] = sizes_pos;
0105       totalSizes_pos += sizes_pos;
0106
0107       int sizes_neg = 0;
0108       for (auto const& mCM_pLS : pLStoLayer[2]) {
0109         std::vector<unsigned int> connectedModuleDetIds_pLS_neg = mCM_pLS.getConnectedModuleDetIds(isuperbin);
0110         connectedModuleDetIds_neg.insert(connectedModuleDetIds_neg.end(),
0111                                          connectedModuleDetIds_pLS_neg.begin(),
0112                                          connectedModuleDetIds_pLS_neg.end());
0113         sizes_neg += connectedModuleDetIds_pLS_neg.size();
0114       }
0115       pixelMapping.connectedPixelsIndexNeg[isuperbin] = totalSizes_neg;
0116       pixelMapping.connectedPixelsSizesNeg[isuperbin] = sizes_neg;
0117       totalSizes_neg += sizes_neg;
0118     }
0119
0120     nPixels = totalSizes + totalSizes_pos + totalSizes_neg;
0121
0122     return {totalSizes,
0123             connectedModuleDetIds,
0124             totalSizes_pos,
0125             connectedModuleDetIds_pos,
0126             totalSizes_neg,
0127             connectedModuleDetIds_neg};
0128   }
0129
0130   inline void fillConnectedModuleArrayExplicit(Modules modules,
0131                                                ModuleMetaData const& mmd,
0132                                                ModuleConnectionMap const& moduleConnectionMap) {
0133     Params_Modules::ArrayU16xMaxConnected* moduleMap = modules.moduleMap();
0134     uint16_t* nConnectedModules = modules.nConnectedModules();
0135
0136     for (auto it = mmd.detIdToIndex.begin(); it != mmd.detIdToIndex.end(); ++it) {
0137       unsigned int detId = it->first;
0138       uint16_t index = it->second;
0139       auto& connectedModules = moduleConnectionMap.getConnectedModuleDetIds(detId);
0140       nConnectedModules[index] = connectedModules.size();
0141       for (uint16_t i = 0; i < nConnectedModules[index]; i++) {
0142         moduleMap[index][i] = mmd.detIdToIndex.at(connectedModules[i]);
0143       }
0144     }
0145   }
0146
0147   inline void fillMapArraysExplicit(Modules modules, ModuleMetaData const& mmd) {
0148     uint16_t* mapIdx = modules.mapIdx();
0149     unsigned int* mapdetId = modules.mapdetId();
0150
0151     unsigned int counter = 0;
0152     for (auto it = mmd.detIdToIndex.begin(); it != mmd.detIdToIndex.end(); ++it) {
0153       unsigned int detId = it->first;
0154       unsigned int index = it->second;
0155       mapIdx[counter] = index;
0156       mapdetId[counter] = detId;
0157       counter++;
0158     }
0159   }
0160
0161   inline void setDerivedQuantities(unsigned int detId,
0162                                    unsigned short& layer,
0163                                    unsigned short& ring,
0164                                    unsigned short& rod,
0165                                    unsigned short& module,
0166                                    unsigned short& subdet,
0167                                    unsigned short& side,
0168                                    float m_x,
0169                                    float m_y,
0170                                    float m_z,
0171                                    float& eta,
0172                                    float& r) {
0173     subdet = (detId & (7 << 25)) >> 25;
0174     side = (subdet == Endcap) ? (detId & (3 << 23)) >> 23 : (detId & (3 << 18)) >> 18;
0175     layer = (subdet == Endcap) ? (detId & (7 << 18)) >> 18 : (detId & (7 << 20)) >> 20;
0176     ring = (subdet == Endcap) ? (detId & (15 << 12)) >> 12 : 0;
0177     module = (detId & (127 << 2)) >> 2;
0178     rod = (subdet == Endcap) ? 0 : (detId & (127 << 10)) >> 10;
0179
0180     r = std::sqrt(m_x * m_x + m_y * m_y + m_z * m_z);
0181     eta = ((m_z > 0) - (m_z < 0)) * std::acosh(r / std::sqrt(m_x * m_x + m_y * m_y));
0182   }
0183
0184   inline void loadCentroidsFromFile(const char* filePath, ModuleMetaData& mmd, uint16_t& nModules) {
0185     std::ifstream ifile(filePath, std::ios::binary);
0186     if (!ifile.is_open()) {
0187       throw std::runtime_error("Unable to open file: " + std::string(filePath));
0188     }
0189
0190     uint16_t counter = 0;
0191     while (!ifile.eof()) {
0192       unsigned int temp_detId;
0193       float module_x, module_y, module_z;
0194       int module_type;
0195
0196       ifile.read(reinterpret_cast<char*>(&temp_detId), sizeof(temp_detId));
0197       ifile.read(reinterpret_cast<char*>(&module_x), sizeof(module_x));
0198       ifile.read(reinterpret_cast<char*>(&module_y), sizeof(module_y));
0199       ifile.read(reinterpret_cast<char*>(&module_z), sizeof(module_z));
0200       ifile.read(reinterpret_cast<char*>(&module_type), sizeof(module_type));
0201
0202       if (ifile) {
0203         mmd.detIdToIndex[temp_detId] = counter;
0204         mmd.module_x[temp_detId] = module_x;
0205         mmd.module_y[temp_detId] = module_y;
0206         mmd.module_z[temp_detId] = module_z;
0207         mmd.module_type[temp_detId] = module_type;
0208         counter++;
0209       } else {
0210         if (!ifile.eof()) {
0211           throw std::runtime_error("Failed to read data for detId: " + std::to_string(temp_detId));
0212         }
0213       }
0214     }
0215
0216     mmd.detIdToIndex[1] = counter;  //pixel module is the last module in the module list
0217     counter++;
0218     nModules = counter;
0219   }
0220
0221   inline std::shared_ptr<ModulesHostCollection> loadModulesFromFile(MapPLStoLayer const& pLStoLayer,
0222                                                                     const char* moduleMetaDataFilePath,
0223                                                                     uint16_t& nModules,
0224                                                                     uint16_t& nLowerModules,
0225                                                                     unsigned int& nPixels,
0226                                                                     PixelMap& pixelMapping,
0227                                                                     const EndcapGeometry& endcapGeometry,
0228                                                                     const TiltedGeometry& tiltedGeometry,
0229                                                                     const ModuleConnectionMap& moduleConnectionMap) {
0230     ModuleMetaData mmd;
0231
0232     loadCentroidsFromFile(moduleMetaDataFilePath, mmd, nModules);
0233
0234     // TODO: this whole section could use some refactoring
0235     auto [totalSizes,
0236           connectedModuleDetIds,
0237           totalSizes_pos,
0238           connectedModuleDetIds_pos,
0239           totalSizes_neg,
0240           connectedModuleDetIds_neg] = getConnectedPixels(nModules, nPixels, pixelMapping, pLStoLayer);
0241
0242     std::array<int, 2> const modules_sizes{{static_cast<int>(nModules), static_cast<int>(nPixels)}};
0243
0244     auto modulesHC = std::make_shared<ModulesHostCollection>(modules_sizes, cms::alpakatools::host());
0245
0246     auto modules_view = modulesHC->view<ModulesSoA>();
0247
0248     // Getting the underlying data pointers
0249     unsigned int* host_detIds = modules_view.detIds();
0250     short* host_layers = modules_view.layers();
0251     short* host_rings = modules_view.rings();
0252     short* host_rods = modules_view.rods();
0253     short* host_modules = modules_view.modules();
0254     short* host_subdets = modules_view.subdets();
0255     short* host_sides = modules_view.sides();
0256     float* host_eta = modules_view.eta();
0257     float* host_r = modules_view.r();
0258     bool* host_isInverted = modules_view.isInverted();
0259     bool* host_isLower = modules_view.isLower();
0260     bool* host_isAnchor = modules_view.isAnchor();
0261     ModuleType* host_moduleType = modules_view.moduleType();
0262     ModuleLayerType* host_moduleLayerType = modules_view.moduleLayerType();
0263     float* host_dxdys = modules_view.dxdys();
0264     float* host_drdzs = modules_view.drdzs();
0265     uint16_t* host_nModules = &modules_view.nModules();
0266     uint16_t* host_nLowerModules = &modules_view.nLowerModules();
0267     uint16_t* host_partnerModuleIndices = modules_view.partnerModuleIndices();
0268     int* host_lstLayers = modules_view.lstLayers();
0269
0270     //reassign detIdToIndex indices here
0271     nLowerModules = (nModules - 1) / 2;
0272     uint16_t lowerModuleCounter = 0;
0273     uint16_t upperModuleCounter = nLowerModules + 1;
0274     //0 to nLowerModules - 1 => only lower modules, nLowerModules - pixel module, nLowerModules + 1 to nModules => upper modules
0275     for (auto it = mmd.detIdToIndex.begin(); it != mmd.detIdToIndex.end(); it++) {
0276       unsigned int detId = it->first;
0277       float m_x = mmd.module_x[detId];
0278       float m_y = mmd.module_y[detId];
0279       float m_z = mmd.module_z[detId];
0280       unsigned int m_t = mmd.module_type[detId];
0281
0282       float eta, r;
0283
0284       uint16_t index;
0285       unsigned short layer, ring, rod, module, subdet, side;
0286       bool isInverted, isLower;
0287       if (detId == 1) {
0288         layer = 0;
0289         ring = 0;
0290         rod = 0;
0291         module = 0;
0292         subdet = 0;
0293         side = 0;
0294         isInverted = false;
0295         isLower = false;
0296         eta = 0;
0297         r = 0;
0298       } else {
0299         setDerivedQuantities(detId, layer, ring, rod, module, subdet, side, m_x, m_y, m_z, eta, r);
0300         isInverted = parseIsInverted(subdet, side, module, layer);
0301         isLower = parseIsLower(isInverted, detId);
0302       }
0303       if (isLower) {
0304         index = lowerModuleCounter;
0305         lowerModuleCounter++;
0306       } else if (detId != 1) {
0307         index = upperModuleCounter;
0308         upperModuleCounter++;
0309       } else {
0310         index = nLowerModules;  //pixel
0311       }
0312       //reassigning indices!
0313       mmd.detIdToIndex[detId] = index;
0314       host_detIds[index] = detId;
0315       host_layers[index] = layer;
0316       host_rings[index] = ring;
0317       host_rods[index] = rod;
0318       host_modules[index] = module;
0319       host_subdets[index] = subdet;
0320       host_sides[index] = side;
0321       host_eta[index] = eta;
0322       host_r[index] = r;
0323       host_isInverted[index] = isInverted;
0324       host_isLower[index] = isLower;
0325
0326       //assigning other variables!
0327       if (detId == 1) {
0328         host_moduleType[index] = PixelModule;
0329         host_moduleLayerType[index] = lst::InnerPixelLayer;
0330         host_dxdys[index] = 0;
0331         host_drdzs[index] = 0;
0332         host_isAnchor[index] = false;
0333       } else {
0334         host_moduleType[index] = (m_t == 25 ? lst::TwoS : lst::PS);
0335         host_moduleLayerType[index] = (m_t == 23 ? lst::Pixel : lst::Strip);
0336
0337         if ((host_moduleType[index] == lst::PS and host_moduleLayerType[index] == lst::Pixel) ||
0338             (host_moduleType[index] == lst::TwoS and host_isLower[index])) {
0339           host_isAnchor[index] = true;
0340         } else {
0341           host_isAnchor[index] = false;
0342         }
0343
0344         host_dxdys[index] = (subdet == Endcap) ? endcapGeometry.getdxdy_slope(detId) : tiltedGeometry.getDxDy(detId);
0345         host_drdzs[index] = (subdet == Barrel) ? tiltedGeometry.getDrDz(detId) : 0;
0346       }
0347
0348       host_lstLayers[index] =
0349           layer + 6 * (subdet == lst::Endcap) + 5 * (subdet == lst::Endcap and host_moduleType[index] == lst::TwoS);
0350     }
0351
0352     //partner module stuff, and slopes and drdz move around
0353     for (auto it = mmd.detIdToIndex.begin(); it != mmd.detIdToIndex.end(); it++) {
0354       auto& detId = it->first;
0355       auto& index = it->second;
0356       if (detId != 1) {
0357         host_partnerModuleIndices[index] =
0358             mmd.detIdToIndex[parsePartnerModuleId(detId, host_isLower[index], host_isInverted[index])];
0359         //add drdz and slope importing stuff here!
0360         if (host_drdzs[index] == 0) {
0361           host_drdzs[index] = host_drdzs[host_partnerModuleIndices[index]];
0362         }
0363         if (host_dxdys[index] == 0) {
0364           host_dxdys[index] = host_dxdys[host_partnerModuleIndices[index]];
0365         }
0366       }
0367     }
0368
0369     *host_nModules = nModules;
0370     *host_nLowerModules = nLowerModules;
0371
0372     // Fill pixel part
0373     pixelMapping.pixelModuleIndex = mmd.detIdToIndex.at(1);
0374
0375     auto modulesPixel_view = modulesHC->view<ModulesPixelSoA>();
0376     auto connectedPixels =
0377         cms::alpakatools::make_host_view(modulesPixel_view.connectedPixels(), modulesPixel_view.metadata().size());
0378     for (unsigned int icondet = 0; icondet < totalSizes; icondet++) {
0379       connectedPixels[icondet] = mmd.detIdToIndex.at(connectedModuleDetIds[icondet]);
0380     }
0381     for (unsigned int icondet = 0; icondet < totalSizes_pos; icondet++) {
0382       connectedPixels[icondet + totalSizes] = mmd.detIdToIndex.at(connectedModuleDetIds_pos[icondet]);
0383     }
0384     for (unsigned int icondet = 0; icondet < totalSizes_neg; icondet++) {
0385       connectedPixels[icondet + totalSizes + totalSizes_pos] = mmd.detIdToIndex.at(connectedModuleDetIds_neg[icondet]);
0386     }
0387
0388     fillConnectedModuleArrayExplicit(modules_view, mmd, moduleConnectionMap);
0389     fillMapArraysExplicit(modules_view, mmd);
0390
0391     return modulesHC;
0392   }
0393 }  // namespace lst
0394 #endif