TrackerTFP/src/LayerEncoding.cc

0001 #include "L1Trigger/TrackerTFP/interface/LayerEncoding.h"
0002
0003 #include <vector>
0004 #include <set>
0005 #include <algorithm>
0006 #include <cmath>
0007 #include <sstream>
0008 #include <fstream>
0009
0010 using namespace std;
0011 using namespace tt;
0012
0013 namespace trackerTFP {
0014
0015   LayerEncoding::LayerEncoding(const DataFormats* dataFormats)
0016       : setup_(dataFormats->setup()),
0017         dataFormats_(dataFormats),
0018         zT_(&dataFormats->format(Variable::zT, Process::zht)),
0019         cot_(&dataFormats->format(Variable::cot, Process::zht)),
0020         layerEncoding_(setup_->numSectorsEta(),
0021                        vector<vector<vector<int>>>(pow(2, zT_->width()), vector<vector<int>>(pow(2, cot_->width())))),
0022         layerEncodingMap_(setup_->numSectorsEta(),
0023                           vector<vector<map<int, const SensorModule*>>>(
0024                               pow(2, zT_->width()), vector<map<int, const SensorModule*>>(pow(2, cot_->width())))),
0025         maybeLayer_(setup_->numSectorsEta(),
0026                     vector<vector<vector<int>>>(pow(2, zT_->width()), vector<vector<int>>(pow(2, cot_->width())))) {
0027     // number of boundaries of fiducial area in r-z plane for a given set of rough r-z track parameter
0028     static constexpr int boundaries = 2;
0029     // find unique sensor mouldes in r-z
0030     // allowed distance in r and z in cm between modules to consider them not unique
0031     static constexpr double delta = 1.e-3;
0032     vector<const SensorModule*> sensorModules;
0033     sensorModules.reserve(setup_->sensorModules().size());
0034     for (const SensorModule& sm : setup_->sensorModules())
0035       sensorModules.push_back(&sm);
0036     auto smallerR = [](const SensorModule* lhs, const SensorModule* rhs) { return lhs->r() < rhs->r(); };
0037     auto smallerZ = [](const SensorModule* lhs, const SensorModule* rhs) { return lhs->z() < rhs->z(); };
0038     auto equalRZ = [](const SensorModule* lhs, const SensorModule* rhs) {
0039       return abs(lhs->r() - rhs->r()) < delta && abs(lhs->z() - rhs->z()) < delta;
0040     };
0041     stable_sort(sensorModules.begin(), sensorModules.end(), smallerR);
0042     stable_sort(sensorModules.begin(), sensorModules.end(), smallerZ);
0043     sensorModules.erase(unique(sensorModules.begin(), sensorModules.end(), equalRZ), sensorModules.end());
0044     // find set of moudles for each set of rough r-z track parameter
0045     // loop over eta sectors
0046     for (int binEta = 0; binEta < setup_->numSectorsEta(); binEta++) {
0047       // cotTheta of eta sector centre
0048       const double sectorCot = (sinh(setup_->boundarieEta(binEta + 1)) + sinh(setup_->boundarieEta(binEta))) / 2.;
0049       // z at radius choenRofZ of eta sector centre
0050       const double sectorZT = setup_->chosenRofZ() * sectorCot;
0051       // loop over bins in zT
0052       for (int binZT = 0; binZT < pow(2, zT_->width()); binZT++) {
0053         // z at radius chosenRofZ wrt zT of sectorZT of this bin centre
0054         const double zT = zT_->floating(zT_->toSigned(binZT));
0055         // z at radius chosenRofZ wrt zT of sectorZT of this bin boundaries
0056         const vector<double> zTs = {sectorZT + zT - zT_->base() / 2., sectorZT + zT + zT_->base() / 2.};
0057         // loop over bins in cotTheta
0058         for (int binCot = 0; binCot < pow(2, cot_->width()); binCot++) {
0059           // cotTheta wrt sectorCot of this bin centre
0060           const double cot = cot_->floating(cot_->toSigned(binCot));
0061           // layer ids crossed by left and right rough r-z parameter shape boundaries
0062           vector<set<int>> layers(boundaries);
0063           map<int, const SensorModule*>& layermaps = layerEncodingMap_[binEta][binZT][binCot];
0064           // cotTheta wrt sectorCot of this bin boundaries
0065           const vector<double> cots = {sectorCot + cot - cot_->base() / 2., sectorCot + cot + cot_->base() / 2.};
0066           // loop over all unique modules
0067           for (const SensorModule* sm : sensorModules) {
0068             // check if module is crossed by left and right rough r-z parameter shape boundaries
0069             for (int i = 0; i < boundaries; i++) {
0070               const int j = boundaries - i - 1;
0071               const double zTi = zTs[sm->r() > setup_->chosenRofZ() ? i : j];
0072               const double coti = cots[sm->r() > setup_->chosenRofZ() ? j : i];
0073               // distance between module and boundary in moudle tilt angle direction
0074               const double d =
0075                   (zTi - sm->z() + (sm->r() - setup_->chosenRofZ()) * coti) / (sm->cosTilt() - sm->sinTilt() * coti);
0076               // compare distance with module size and add module layer id to layers if module is crossed
0077               if (abs(d) < sm->numColumns() * sm->pitchCol() / 2.) {
0078                 layers[i].insert(sm->layerId());
0079                 layermaps[sm->layerId()] = sm;
0080               }
0081             }
0082           }
0083           // mayber layers are given by layer ids crossed by only one booundary
0084           set<int> maybeLayer;
0085           set_symmetric_difference(layers[0].begin(),
0086                                    layers[0].end(),
0087                                    layers[1].begin(),
0088                                    layers[1].end(),
0089                                    inserter(maybeLayer, maybeLayer.end()));
0090           // layerEncoding is given by sorted layer ids crossed by any booundary
0091           set<int> layerEncoding;
0092           set_union(layers[0].begin(),
0093                     layers[0].end(),
0094                     layers[1].begin(),
0095                     layers[1].end(),
0096                     inserter(layerEncoding, layerEncoding.end()));
0097           vector<int>& le = layerEncoding_[binEta][binZT][binCot];
0098           le = vector<int>(layerEncoding.begin(), layerEncoding.end());
0099           vector<int>& ml = maybeLayer_[binEta][binZT][binCot];
0100           ml.reserve(maybeLayer.size());
0101           for (int m : maybeLayer) {
0102             int layer = distance(le.begin(), find(le.begin(), le.end(), m));
0103             if (layer >= setup_->numLayers())
0104               layer = setup_->numLayers() - 1;
0105             ml.push_back(layer);
0106           }
0107         }
0108       }
0109     }
0110     const bool print = false;
0111     if (!print)
0112       return;
0113     static constexpr int widthLayer = 3;
0114     static constexpr auto layerIds = {1, 2, 3, 4, 5, 6, 11, 12, 13, 14, 15};
0115     stringstream ss;
0116     for (int layer : layerIds) {
0117       auto encode = [layer, this](const vector<int>& layers, int& l) {
0118         const auto it = find(layers.begin(), layers.end(), layer);
0119         if (it == layers.end())
0120           return false;
0121         l = distance(layers.begin(), it);
0122         if (l >= setup_->numLayers())
0123           l = setup_->numLayers() - 1;
0124         return true;
0125       };
0126       for (int binEta = 0; binEta < setup_->numSectorsEta(); binEta++) {
0127         for (int binZT = 0; binZT < pow(2, zT_->width()); binZT++) {
0128           for (int binCot = 0; binCot < pow(2, cot_->width()); binCot++) {
0129             const int zT =
0130                 binZT < pow(2, zT_->width() - 1) ? binZT + pow(2, zT_->width() - 1) : binZT - pow(2, zT_->width() - 1);
0131             const int cot = binCot < pow(2, cot_->width() - 1) ? binCot + pow(2, cot_->width() - 1)
0132                                                                : binCot - pow(2, cot_->width() - 1);
0133             const vector<int>& layers = layerEncoding_[binEta][zT][cot];
0134             const vector<int>& maybes = maybeLayer_[binEta][zT][cot];
0135             int layerKF(-1);
0136             if (encode(layers, layerKF))
0137               ss << "1" << TTBV(layerKF, widthLayer) << (encode(maybes, layerKF) ? "1" : "0");
0138             else
0139               ss << "00000";
0140             ss << endl;
0141           }
0142         }
0143       }
0144     }
0145     fstream file;
0146     file.open("layerEncoding.txt", ios::out);
0147     file << ss.rdbuf();
0148     file.close();
0149   }
0150
0151   // encoded layer id for given eta sector, bin in zT, bin in cotThea and decoed layer id, returns -1 if layer incositent with track
0152   const int LayerEncoding::layerIdKF(int binEta, int binZT, int binCot, int layerId) const {
0153     const vector<int>& layers = layerEncoding_[binEta][binZT][binCot];
0154     const auto it = find(layers.begin(), layers.end(), layerId);
0155     if (it == layers.end())
0156       return -1;
0157     int layer = distance(layers.begin(), it);
0158     if (layer >= setup_->numLayers())
0159       layer = setup_->numLayers() - 1;
0160     return layer;
0161   }
0162
0163   // pattern of maybe layers for given eta sector, bin in zT and bin in cotThea
0164   TTBV LayerEncoding::maybePattern(int binEta, int binZT, int binCot) const {
0165     TTBV ttBV(0, setup_->numLayers());
0166     const vector<int>& layers = layerEncoding_[binEta][binZT][binCot];
0167     const vector<int>& maybes = maybeLayer_[binEta][binZT][binCot];
0168     for (int m : maybes)
0169       ttBV.set(distance(layers.begin(), find(layers.begin(), layers.end(), m)));
0170     return ttBV;
0171   }
0172
0173 }  // namespace trackerTFP