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 namespace trackerTFP {
0011
0012   LayerEncoding::LayerEncoding(const DataFormats* dataFormats)
0013       : setup_(dataFormats->setup()),
0014         dataFormats_(dataFormats),
0015         zT_(&dataFormats->format(Variable::zT, Process::gp)),
0016         layerEncoding_(std::vector<std::vector<int>>(pow(2, zT_->width()))),
0017         maybePattern_(std::vector<TTBV>(pow(2, zT_->width()), TTBV(0, setup_->numLayers()))),
0018         nullLE_(setup_->numLayers(), 0),
0019         nullMP_(0, setup_->numLayers()) {
0020     // number of boundaries of fiducial area in r-z plane for a given set of rough r-z track parameter
0021     static constexpr int boundaries = 2;
0022     // z at radius chosenRofZ wrt zT of sectorZT of this bin boundaries
0023     const std::vector<double> z0s = {-setup_->beamWindowZ(), setup_->beamWindowZ()};
0024     // find unique sensor mouldes in r-z
0025     // allowed distance in r and z in cm between modules to consider them not unique
0026     static constexpr double delta = 1.e-3;
0027     std::vector<const tt::SensorModule*> sensorModules;
0028     sensorModules.reserve(setup_->sensorModules().size());
0029     for (const tt::SensorModule& sm : setup_->sensorModules())
0030       sensorModules.push_back(&sm);
0031     auto smallerR = [](const tt::SensorModule* lhs, const tt::SensorModule* rhs) { return lhs->r() < rhs->r(); };
0032     auto smallerZ = [](const tt::SensorModule* lhs, const tt::SensorModule* rhs) { return lhs->z() < rhs->z(); };
0033     auto equalRZ = [](const tt::SensorModule* lhs, const tt::SensorModule* rhs) {
0034       return std::abs(lhs->r() - rhs->r()) < delta && std::abs(lhs->z() - rhs->z()) < delta;
0035     };
0036     std::stable_sort(sensorModules.begin(), sensorModules.end(), smallerZ);
0037     std::stable_sort(sensorModules.begin(), sensorModules.end(), smallerR);
0038     sensorModules.erase(std::unique(sensorModules.begin(), sensorModules.end(), equalRZ), sensorModules.end());
0039     std::stable_sort(sensorModules.begin(), sensorModules.end(), smallerZ);
0040     sensorModules.erase(std::unique(sensorModules.begin(), sensorModules.end(), equalRZ), sensorModules.end());
0041     // find set of moudles for each set of rough r-z track parameter
0042     // loop over zT bins
0043     for (int binZT = 0; binZT < std::pow(2, zT_->width()); binZT++) {
0044       // z at radius chosenRofZ
0045       const double zT = zT_->floating(zT_->toSigned(binZT));
0046       // z at radius chosenRofZ wrt zT of sectorZT of this bin boundaries
0047       const std::vector<double> zTs = {zT - zT_->base() / 2., zT + zT_->base() / 2.};
0048       std::vector<std::vector<double>> cots(boundaries);
0049       for (int i = 0; i < boundaries; i++)
0050         for (double z0 : z0s)
0051           cots[i].push_back((zTs[i] - z0) / setup_->chosenRofZ());
0052       // layer ids crossed by left and right rough r-z parameter shape boundaries
0053       std::vector<std::set<int>> layers(boundaries);
0054       // loop over all unique modules
0055       for (const tt::SensorModule* sm : sensorModules) {
0056         // check if module is crossed by left and right rough r-z parameter shape boundaries
0057         for (int i = 0; i < boundaries; i++) {
0058           const double zTi = zTs[i];
0059           const double coti = sm->r() < setup_->chosenRofZ() ? cots[i][i == 0 ? 0 : 1] : cots[i][i == 0 ? 1 : 0];
0060           // distance between module and boundary in moudle tilt angle direction
0061           const double d =
0062               (zTi - sm->z() + (sm->r() - setup_->chosenRofZ()) * coti) / (sm->cosTilt() - sm->sinTilt() * coti);
0063           // compare distance with module size and add module layer id to layers if module is crossed
0064           if (std::abs(d) < sm->numColumns() * sm->pitchCol() / 2.)
0065             layers[i].insert(sm->layerId());
0066         }
0067       }
0068       // mayber layers are given by layer ids crossed by only one boundary
0069       std::set<int> maybeLayer;
0070       std::set_symmetric_difference(layers[0].begin(),
0071                                     layers[0].end(),
0072                                     layers[1].begin(),
0073                                     layers[1].end(),
0074                                     std::inserter(maybeLayer, maybeLayer.end()));
0075       // layerEncoding is given by sorted layer ids crossed by any boundary
0076       std::set<int> layerEncoding;
0077       std::set_union(layers[0].begin(),
0078                      layers[0].end(),
0079                      layers[1].begin(),
0080                      layers[1].end(),
0081                      std::inserter(layerEncoding, layerEncoding.end()));
0082       // fill layerEncoding_
0083       std::vector<int>& le = layerEncoding_[binZT];
0084       le = std::vector<int>(layerEncoding.begin(), layerEncoding.end());
0085       le.resize(setup_->numLayers(), -1);
0086       // fill maybePattern_
0087       TTBV& mp = maybePattern_[binZT];
0088       for (int m : maybeLayer)
0089         mp.set(std::min(static_cast<int>(std::distance(le.begin(), std::find(le.begin(), le.end(), m))),
0090                         setup_->numLayers() - 1));
0091     }
0092   }
0093
0094   // Set of layers for given bin in zT
0095   const std::vector<int>& LayerEncoding::layerEncoding(int zT) const {
0096     const int binZT = zT_->toUnsigned(zT);
0097     return zT_->inRange(zT) ? layerEncoding_.at(binZT) : nullLE_;
0098   }
0099
0100   // Set of layers for given zT in cm
0101   const std::vector<int>& LayerEncoding::layerEncoding(double zT) const {
0102     const int binZT = zT_->integer(zT);
0103     return layerEncoding(binZT);
0104   }
0105
0106   // pattern of maybe layers for given bin in zT
0107   const TTBV& LayerEncoding::maybePattern(int zT) const {
0108     const int binZT = zT_->toUnsigned(zT);
0109     return zT_->inRange(zT) ? maybePattern_[binZT] : nullMP_;
0110   }
0111
0112   // pattern of maybe layers for given zT in cm
0113   const TTBV& LayerEncoding::maybePattern(double zT) const {
0114     const int binZT = zT_->integer(zT);
0115     return maybePattern(binZT);
0116   }
0117
0118   // fills numPS, num2S, numMissingPS and numMissingPS for given hitPattern and trajectory
0119   void LayerEncoding::analyze(
0120       int hitpattern, double cot, double z0, int& numPS, int& num2S, int& numMissingPS, int& numMissing2S) const {
0121     // look up layer encoding nad maybe pattern
0122     const double zT = z0 + setup_->chosenRofZ() * cot;
0123     const std::vector<int>& le = this->layerEncoding(zT);
0124     const TTBV& mp = this->maybePattern(zT);
0125     const TTBV hp(hitpattern, setup_->numLayers());
0126     // loop from innermost layer to outermost hitted layer
0127     for (int layerIdKF = 0; layerIdKF <= hp.pmEncode(); layerIdKF++) {
0128       // look up layer Id [1-6 barrel, 11-15 disks]
0129       const int layerId = le[layerIdKF];
0130       // identify module type
0131       bool ps = layerId <= setup_->numBarrelLayerPS();
0132       const bool barrel = layerId <= setup_->numBarrelLayer();
0133       if (!barrel) {
0134         // calc disk id (0 - 4)
0135         const int diskId = layerId - setup_->offsetLayerDisks() - setup_->offsetLayerId();
0136         // avergae disk z position
0137         const double z = setup_->hybridDiskZ(diskId) * (cot < 0. ? -1. : 1.);
0138         // innermost edge of 2S modules
0139         const double rLimit = setup_->disk2SR(diskId, 0) - setup_->pitchCol2S();
0140         // trajectory radius at avergae disk z position
0141         const double r = (z - z0) / cot;
0142         // compare with innermost edge of 2S modules to identify PS
0143         if (r < rLimit)
0144           ps = true;
0145       }
0146       if (hp.test(layerIdKF))  // layer is hit
0147         ps ? numPS++ : num2S++;
0148       else if (!mp.test(layerIdKF))  // layer is not hit but should have been hitted (roughly by) trajectory
0149         ps ? numMissingPS++ : numMissing2S++;
0150     }
0151   }
0152
0153 }  // namespace trackerTFP