TrackerTFP/src/LayerEncoding.cc

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