TrackerDTC/src/Stub.cc

0001 #include "L1Trigger/TrackerDTC/interface/Stub.h"
0002
0003 #include <cmath>
0004 #include <iterator>
0005 #include <algorithm>
0006 #include <utility>
0007
0008 using namespace edm;
0009 using namespace std;
0010 using namespace tt;
0011
0012 namespace trackerDTC {
0013
0014   Stub::Stub(const ParameterSet& iConfig,
0015              const Setup* setup,
0016              const LayerEncoding* layerEncoding,
0017              SensorModule* sm,
0018              const TTStubRef& ttStubRef)
0019       : setup_(setup),
0020         layerEncoding_(layerEncoding),
0021         sm_(sm),
0022         ttStubRef_(ttStubRef),
0023         hybrid_(iConfig.getParameter<bool>("UseHybrid")),
0024         valid_(true) {
0025     regions_.reserve(setup->numOverlappingRegions());
0026     // get stub local coordinates
0027     const MeasurementPoint& mp = ttStubRef->clusterRef(0)->findAverageLocalCoordinatesCentered();
0028
0029     // convert to uniformed local coordinates
0030
0031     // column number in pitch units
0032     col_ = (int)floor(pow(-1, sm->signCol()) * (mp.y() - sm->numColumns() / 2) / setup->baseCol());
0033     // row number in half pitch units
0034     row_ = (int)floor(pow(-1, sm->signRow()) * (mp.x() - sm->numRows() / 2) / setup->baseRow());
0035     // bend number in quarter pitch units
0036     bend_ = (int)floor(pow(-1, sm->signBend()) * (ttStubRef->bendBE()) / setup->baseBend());
0037     // reduced row number for look up
0038     rowLUT_ = (int)floor((double)row_ / pow(2., setup->widthRow() - setup->dtcWidthRowLUT()));
0039     // sub row number inside reduced row number
0040     rowSub_ = row_ - (rowLUT_ + .5) * pow(2, setup->widthRow() - setup->dtcWidthRowLUT());
0041
0042     // convert local to global coordinates
0043
0044     const double y = (col_ + .5) * setup->baseCol() * sm->pitchCol();
0045     // radius of a column of strips/pixel in cm
0046     d_ = sm->r() + y * sm->sinTilt();
0047     // stub z in cm
0048     z_ = digi(sm->z() + y * sm->cosTilt(), setup->tmttBaseZ());
0049
0050     const double x0 = rowLUT_ * setup->baseRow() * setup->dtcNumMergedRows() * sm->pitchRow();
0051     const double x1 = (rowLUT_ + 1) * setup->baseRow() * setup->dtcNumMergedRows() * sm->pitchRow();
0052     const double x = (rowLUT_ + .5) * setup->baseRow() * setup->dtcNumMergedRows() * sm->pitchRow();
0053     // stub r in cm
0054     r_ = sqrt(d_ * d_ + x * x);
0055
0056     const double phi0 = sm->phi() + atan2(x0, d_);
0057     const double phi1 = sm->phi() + atan2(x1, d_);
0058     const double c = (phi0 + phi1) / 2.;
0059     const double m = (phi1 - phi0) / setup->dtcNumMergedRows();
0060
0061     // intercept of linearized stub phi in rad
0062     c_ = digi(c, setup->tmttBasePhi());
0063     // slope of linearized stub phi in rad / strip
0064     m_ = digi(m, setup->dtcBaseM());
0065
0066     if (hybrid_) {
0067       if (abs(z_ / r_) > setup->hybridMaxCot())
0068         // did not pass eta cut
0069         valid_ = false;
0070     } else {
0071       // extrapolated z at radius T assuming z0=0
0072       const double zT = setup->chosenRofZ() * z_ / r_;
0073       // extrapolated z0 window at radius T
0074       const double dZT = setup->beamWindowZ() * abs(1. - setup->chosenRofZ() / r_);
0075       double zTMin = zT - dZT;
0076       double zTMax = zT + dZT;
0077       if (zTMin >= setup->maxZT() || zTMax < -setup->maxZT())
0078         // did not pass "eta" cut
0079         valid_ = false;
0080       else {
0081         zTMin = max(zTMin, -setup->maxZT());
0082         zTMax = min(zTMax, setup->maxZT());
0083       }
0084       // range of stub cot(theta)
0085       cot_ = {zTMin / setup->chosenRofZ(), zTMax / setup->chosenRofZ()};
0086     }
0087
0088     // stub r w.r.t. chosenRofPhi in cm
0089     static const double chosenRofPhi = hybrid_ ? setup->hybridChosenRofPhi() : setup->chosenRofPhi();
0090     r_ = digi(r_ - chosenRofPhi, setup->tmttBaseR());
0091
0092     // radial (cylindrical) component of sensor separation
0093     const double dr = sm->sep() / (sm->cosTilt() - sm->sinTilt() * z_ / d_);
0094     // converts bend into inv2R in 1/cm
0095     const double inv2ROverBend = sm->pitchRow() / dr / d_;
0096     // inv2R in 1/cm
0097     const double inv2R = -bend_ * setup->baseBend() * inv2ROverBend;
0098     // inv2R uncertainty in 1/cm
0099     const double dInv2R = setup->bendCut() * inv2ROverBend;
0100     const double minPt = hybrid_ ? setup->hybridMinPtStub() : setup->minPt();
0101     const double maxInv2R = setup->invPtToDphi() / minPt - setup->dtcBaseInv2R() / 2.;
0102     double inv2RMin = digi(inv2R - dInv2R, setup->dtcBaseInv2R());
0103     double inv2RMax = digi(inv2R + dInv2R, setup->dtcBaseInv2R());
0104     if (inv2RMin > maxInv2R || inv2RMax < -maxInv2R) {
0105       // did not pass pt cut
0106       valid_ = false;
0107     } else {
0108       inv2RMin = max(inv2RMin, -maxInv2R);
0109       inv2RMax = min(inv2RMax, maxInv2R);
0110     }
0111     // range of stub inv2R in 1/cm
0112     inv2R_ = {inv2RMin, inv2RMax};
0113
0114     // stub phi w.r.t. detector region centre in rad
0115     phi_ = c_ + rowSub_ * m_;
0116
0117     // range of stub extrapolated phi to radius chosenRofPhi in rad
0118     phiT_.first = phi_ - r_ * inv2R_.first;
0119     phiT_.second = phi_ - r_ * inv2R_.second;
0120     if (phiT_.first > phiT_.second)
0121       swap(phiT_.first, phiT_.second);
0122
0123     if (phiT_.first < 0.)
0124       regions_.push_back(0);
0125     if (phiT_.second >= 0.)
0126       regions_.push_back(1);
0127
0128     // apply data format specific manipulations
0129     if (!hybrid_)
0130       return;
0131
0132     // stub r w.r.t. an offset in cm
0133     r_ -= sm->offsetR() - chosenRofPhi;
0134     // stub z w.r.t. an offset in cm
0135     z_ -= sm->offsetZ();
0136     if (sm->type() == SensorModule::Disk2S) {
0137       // encoded r
0138       r_ = sm->encodedR() + (sm->side() ? -col_ : (col_ + sm->numColumns() / 2));
0139       r_ = (r_ + 0.5) * setup->hybridBaseR(sm->type());
0140     }
0141
0142     // encode bend
0143     const vector<double>& encodingBend = setup->encodingBend(sm->windowSize(), sm->psModule());
0144     const auto pos = find(encodingBend.begin(), encodingBend.end(), abs(ttStubRef->bendBE()));
0145     const int uBend = distance(encodingBend.begin(), pos);
0146     bend_ = pow(-1, signbit(bend_)) * uBend;
0147   }
0148
0149   // returns bit accurate representation of Stub
0150   Frame Stub::frame(int region) const { return hybrid_ ? formatHybrid(region) : formatTMTT(region); }
0151
0152   // returns true if stub belongs to region
0153   bool Stub::inRegion(int region) const { return find(regions_.begin(), regions_.end(), region) != regions_.end(); }
0154
0155   // truncates double precision to f/w integer equivalent
0156   double Stub::digi(double value, double precision) const { return (floor(value / precision) + .5) * precision; }
0157
0158   // returns 64 bit stub in hybrid data format
0159   Frame Stub::formatHybrid(int region) const {
0160     const SensorModule::Type type = sm_->type();
0161     // layer encoding
0162     const int decodedLayerId = layerEncoding_->decode(sm_);
0163     // stub phi w.r.t. processing region border in rad
0164     double phi = phi_ - (region - .5) * setup_->baseRegion() + setup_->hybridRangePhi() / 2.;
0165     if (phi >= setup_->hybridRangePhi())
0166       phi = setup_->hybridRangePhi() - setup_->hybridBasePhi(type) / 2.;
0167     // convert stub variables into bit vectors
0168     const TTBV hwR(r_, setup_->hybridBaseR(type), setup_->hybridWidthR(type), true);
0169     const TTBV hwPhi(phi, setup_->hybridBasePhi(type), setup_->hybridWidthPhi(type), true);
0170     const TTBV hwZ(z_, setup_->hybridBaseZ(type), setup_->hybridWidthZ(type), true);
0171     const TTBV hwAlpha(row_, setup_->hybridBaseAlpha(type), setup_->hybridWidthAlpha(type), true);
0172     const TTBV hwBend(bend_, setup_->hybridWidthBend(type), true);
0173     const TTBV hwLayer(decodedLayerId, setup_->hybridWidthLayerId());
0174     const TTBV hwGap(0, setup_->hybridNumUnusedBits(type));
0175     const TTBV hwValid(1, 1);
0176     // assemble final bitset
0177     return Frame(hwGap.str() + hwR.str() + hwZ.str() + hwPhi.str() + hwAlpha.str() + hwBend.str() + hwLayer.str() +
0178                  hwValid.str());
0179   }
0180
0181   Frame Stub::formatTMTT(int region) const {
0182     int layerM = sm_->layerId();
0183     // convert unique layer id [1-6,11-15] into reduced layer id [0-6]
0184     // a fiducial track may not cross more then 7 detector layers, for stubs from a given track the reduced layer id is actually unique
0185     int layer(-1);
0186     if (layerM == 1)
0187       layer = 0;
0188     else if (layerM == 2)
0189       layer = 1;
0190     else if (layerM == 6 || layerM == 11)
0191       layer = 2;
0192     else if (layerM == 5 || layerM == 12)
0193       layer = 3;
0194     else if (layerM == 4 || layerM == 13)
0195       layer = 4;
0196     else if (layerM == 14)
0197       layer = 5;
0198     else if (layerM == 3 || layerM == 15)
0199       layer = 6;
0200     // assign stub to phi sectors within a processing region, to be generalized
0201     TTBV sectorsPhi(0, setup_->numOverlappingRegions() * setup_->numSectorsPhi());
0202     if (phiT_.first < 0.) {
0203       if (phiT_.first < -setup_->baseSector())
0204         sectorsPhi.set(0);
0205       else
0206         sectorsPhi.set(1);
0207       if (phiT_.second < 0. && phiT_.second >= -setup_->baseSector())
0208         sectorsPhi.set(1);
0209     }
0210     if (phiT_.second >= 0.) {
0211       if (phiT_.second < setup_->baseSector())
0212         sectorsPhi.set(2);
0213       else
0214         sectorsPhi.set(3);
0215       if (phiT_.first >= 0. && phiT_.first < setup_->baseSector())
0216         sectorsPhi.set(2);
0217     }
0218     // assign stub to eta sectors within a processing region
0219     pair<int, int> sectorEta({0, setup_->numSectorsEta() - 1});
0220     for (int bin = 0; bin < setup_->numSectorsEta(); bin++)
0221       if (asinh(cot_.first) < setup_->boundarieEta(bin + 1)) {
0222         sectorEta.first = bin;
0223         break;
0224       }
0225     for (int bin = sectorEta.first; bin < setup_->numSectorsEta(); bin++)
0226       if (asinh(cot_.second) < setup_->boundarieEta(bin + 1)) {
0227         sectorEta.second = bin;
0228         break;
0229       }
0230     // stub phi w.r.t. processing region centre in rad
0231     const double phi = phi_ - (region - .5) * setup_->baseRegion();
0232     // convert stub variables into bit vectors
0233     const TTBV hwValid(1, 1);
0234     const TTBV hwGap(0, setup_->tmttNumUnusedBits());
0235     const TTBV hwLayer(layer, setup_->tmttWidthLayer());
0236     const TTBV hwSectorEtaMin(sectorEta.first, setup_->tmttWidthSectorEta());
0237     const TTBV hwSectorEtaMax(sectorEta.second, setup_->tmttWidthSectorEta());
0238     const TTBV hwR(r_, setup_->tmttBaseR(), setup_->tmttWidthR(), true);
0239     const TTBV hwPhi(phi, setup_->tmttBasePhi(), setup_->tmttWidthPhi(), true);
0240     const TTBV hwZ(z_, setup_->tmttBaseZ(), setup_->tmttWidthZ(), true);
0241     const TTBV hwInv2RMin(inv2R_.first, setup_->tmttBaseInv2R(), setup_->tmttWidthInv2R(), true);
0242     const TTBV hwInv2RMax(inv2R_.second, setup_->tmttBaseInv2R(), setup_->tmttWidthInv2R(), true);
0243     TTBV hwSectorPhis(0, setup_->numSectorsPhi());
0244     for (int sectorPhi = 0; sectorPhi < setup_->numSectorsPhi(); sectorPhi++)
0245       hwSectorPhis[sectorPhi] = sectorsPhi[region * setup_->numSectorsPhi() + sectorPhi];
0246     // assemble final bitset
0247     return Frame(hwGap.str() + hwValid.str() + hwR.str() + hwPhi.str() + hwZ.str() + hwLayer.str() +
0248                  hwSectorPhis.str() + hwSectorEtaMin.str() + hwSectorEtaMax.str() + hwInv2RMin.str() +
0249                  hwInv2RMax.str());
0250   }
0251
0252 }  // namespace trackerDTC