CSCTrackFinder/src/CSCSectorReceiverLUT.cc

0001 #include "L1Trigger/CSCTrackFinder/interface/CSCSectorReceiverLUT.h"
0002 #include "L1Trigger/CSCTrackFinder/interface/CSCSectorReceiverMiniLUT.h"
0003 #include "L1Trigger/CSCTriggerPrimitives/interface/CSCPatternBank.h"
0004 #include "DataFormats/L1CSCTrackFinder/interface/CSCBitWidths.h"
0005 #include "DataFormats/L1CSCTrackFinder/interface/CSCTFConstants.h"
0006 #include "DataFormats/CSCDigi/interface/CSCConstants.h"
0007
0008 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0009 #include "Geometry/CSCGeometry/interface/CSCLayerGeometry.h"
0010 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0011 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0012
0013 #include "DataFormats/MuonDetId/interface/CSCTriggerNumbering.h"
0014
0015 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0016
0017 #include <fstream>
0018 #include <cstring>
0019
0020 lclphidat* CSCSectorReceiverLUT::me_lcl_phi = nullptr;
0021 bool CSCSectorReceiverLUT::me_lcl_phi_loaded = false;
0022
0023 CSCSectorReceiverLUT::CSCSectorReceiverLUT(
0024     int endcap, int sector, int subsector, int station, const edm::ParameterSet& pset, bool TMB07)
0025     : _endcap(endcap), _sector(sector), _subsector(subsector), _station(station), isTMB07(TMB07) {
0026   LUTsFromFile = pset.getUntrackedParameter<bool>("ReadLUTs", false);
0027   useMiniLUTs = pset.getUntrackedParameter<bool>("UseMiniLUTs", true);
0028   isBinary = pset.getUntrackedParameter<bool>("Binary", false);
0029
0030   me_global_eta = nullptr;
0031   me_global_phi = nullptr;
0032   mb_global_phi = nullptr;
0033   if (LUTsFromFile && !useMiniLUTs) {
0034     me_lcl_phi_file = pset.getUntrackedParameter<edm::FileInPath>(
0035         "LocalPhiLUT",
0036         edm::FileInPath(std::string("L1Trigger/CSCTrackFinder/LUTs/LocalPhiLUT" +
0037                                     (isBinary ? std::string(".bin") : std::string(".dat")))));
0038     me_gbl_phi_file = pset.getUntrackedParameter<edm::FileInPath>(
0039         "GlobalPhiLUTME",
0040         edm::FileInPath((std::string("L1Trigger/CSCTrackFinder/LUTs/GlobalPhiME") + encodeFileIndex() +
0041                          (isBinary ? std::string(".bin") : std::string(".dat")))));
0042     if (station == 1)
0043       mb_gbl_phi_file = pset.getUntrackedParameter<edm::FileInPath>(
0044           "GlobalPhiLUTMB",
0045           edm::FileInPath((std::string("L1Trigger/CSCTrackFinder/LUTs/GlobalPhiMB") + encodeFileIndex() +
0046                            (isBinary ? std::string(".bin") : std::string(".dat")))));
0047     me_gbl_eta_file = pset.getUntrackedParameter<edm::FileInPath>(
0048         "GlobalEtaLUTME",
0049         edm::FileInPath((std::string("L1Trigger/CSCTrackFinder/LUTs/GlobalEtaME") + encodeFileIndex() +
0050                          (isBinary ? std::string(".bin") : std::string(".dat")))));
0051     readLUTsFromFile();
0052   }
0053 }
0054
0055 CSCSectorReceiverLUT::CSCSectorReceiverLUT(const CSCSectorReceiverLUT& lut)
0056     : _endcap(lut._endcap),
0057       _sector(lut._sector),
0058       _subsector(lut._subsector),
0059       _station(lut._station),
0060       me_lcl_phi_file(lut.me_lcl_phi_file),
0061       me_gbl_phi_file(lut.me_gbl_phi_file),
0062       mb_gbl_phi_file(lut.mb_gbl_phi_file),
0063       me_gbl_eta_file(lut.me_gbl_eta_file),
0064       LUTsFromFile(lut.LUTsFromFile),
0065       isBinary(lut.isBinary) {
0066   if (lut.mb_global_phi) {
0067     mb_global_phi = new gblphidat[1 << CSCBitWidths::kGlobalPhiAddressWidth];
0068     memcpy(mb_global_phi, lut.mb_global_phi, (1 << CSCBitWidths::kGlobalPhiAddressWidth) * sizeof(gblphidat));
0069   } else
0070     mb_global_phi = nullptr;
0071   if (lut.me_global_phi) {
0072     me_global_phi = new gblphidat[1 << CSCBitWidths::kGlobalPhiAddressWidth];
0073     memcpy(me_global_phi, lut.me_global_phi, (1 << CSCBitWidths::kGlobalPhiAddressWidth) * sizeof(gblphidat));
0074   } else
0075     me_global_phi = nullptr;
0076   if (lut.me_global_eta) {
0077     me_global_eta = new gbletadat[1 << CSCBitWidths::kGlobalEtaAddressWidth];
0078     memcpy(me_global_eta, lut.me_global_eta, (1 << CSCBitWidths::kGlobalEtaAddressWidth) * sizeof(gbletadat));
0079   } else
0080     me_global_eta = nullptr;
0081 }
0082
0083 CSCSectorReceiverLUT& CSCSectorReceiverLUT::operator=(const CSCSectorReceiverLUT& lut) {
0084   if (this != &lut) {
0085     _endcap = lut._endcap;
0086     _sector = lut._sector;
0087     _subsector = lut._subsector;
0088     _station = lut._station;
0089     me_lcl_phi_file = lut.me_lcl_phi_file;
0090     me_gbl_phi_file = lut.me_gbl_phi_file;
0091     mb_gbl_phi_file = lut.mb_gbl_phi_file;
0092     me_gbl_eta_file = lut.me_gbl_eta_file;
0093     LUTsFromFile = lut.LUTsFromFile;
0094     isBinary = lut.isBinary;
0095
0096     if (lut.mb_global_phi) {
0097       mb_global_phi = new gblphidat[1 << CSCBitWidths::kGlobalPhiAddressWidth];
0098       memcpy(mb_global_phi, lut.mb_global_phi, (1 << CSCBitWidths::kGlobalPhiAddressWidth) * sizeof(gblphidat));
0099     } else
0100       mb_global_phi = nullptr;
0101
0102     if (lut.me_global_phi) {
0103       me_global_phi = new gblphidat[1 << CSCBitWidths::kGlobalPhiAddressWidth];
0104       memcpy(me_global_phi, lut.me_global_phi, (1 << CSCBitWidths::kGlobalPhiAddressWidth) * sizeof(gblphidat));
0105     } else
0106       me_global_phi = nullptr;
0107
0108     if (lut.me_global_eta) {
0109       me_global_eta = new gbletadat[1 << CSCBitWidths::kGlobalEtaAddressWidth];
0110       memcpy(me_global_eta, lut.me_global_eta, (1 << CSCBitWidths::kGlobalEtaAddressWidth) * sizeof(gbletadat));
0111     } else
0112       me_global_eta = nullptr;
0113   }
0114   return *this;
0115 }
0116
0117 CSCSectorReceiverLUT::~CSCSectorReceiverLUT() {
0118   if (me_lcl_phi_loaded) {
0119     delete me_lcl_phi;
0120     me_lcl_phi = nullptr;
0121     me_lcl_phi_loaded = false;
0122   }
0123   if (me_global_eta) {
0124     delete me_global_eta;
0125     me_global_eta = nullptr;
0126   }
0127   if (me_global_phi) {
0128     delete me_global_phi;
0129     me_global_phi = nullptr;
0130   }
0131   if (mb_global_phi) {
0132     delete mb_global_phi;
0133     mb_global_phi = nullptr;
0134   }
0135 }
0136
0137 lclphidat CSCSectorReceiverLUT::calcLocalPhi(const lclphiadd& theadd) const {
0138   lclphidat data;
0139
0140   constexpr int maxPhiL = 1 << CSCBitWidths::kLocalPhiDataBitWidth;
0141   double binPhiL = static_cast<double>(maxPhiL) / (2. * CSCConstants::MAX_NUM_STRIPS_RUN1);
0142
0143   double patternOffset;
0144
0145   patternOffset = CSCPatternBank::getLegacyPosition((theadd.pattern_type << 3) + theadd.clct_pattern);
0146
0147   // The phiL value stored is for the center of the half-/di-strip.
0148   if (theadd.strip < 2 * CSCConstants::MAX_NUM_STRIPS_RUN1)
0149     if (theadd.pattern_type == 1 || isTMB07)  // if halfstrip (Note: no distrips in TMB 2007 patterns)
0150       data.phi_local = static_cast<unsigned>((0.5 + theadd.strip + patternOffset) * binPhiL);
0151     else  // if distrip
0152       data.phi_local = static_cast<unsigned>((2 + theadd.strip + 4. * patternOffset) * binPhiL);
0153   else {
0154     throw cms::Exception("CSCSectorReceiverLUT") << "+++ Value of strip, " << theadd.strip << ", exceeds max allowed, "
0155                                                  << 2 * CSCConstants::MAX_NUM_STRIPS_RUN1 - 1 << " +++\n";
0156   }
0157
0158   if (data.phi_local >= maxPhiL) {
0159     throw cms::Exception("CSCSectorReceiverLUT")
0160         << "+++ Value of phi_local, " << data.phi_local << ", exceeds max allowed, " << maxPhiL - 1 << " +++\n";
0161   }
0162
0163   LogDebug("CSCSectorReceiver") << "endcap = " << _endcap << " station = " << _station << " maxPhiL = " << maxPhiL
0164                                 << " binPhiL = " << binPhiL;
0165   LogDebug("CSCSectorReceiver") << "strip # " << theadd.strip << " hs/ds = " << theadd.pattern_type
0166                                 << " pattern = " << theadd.clct_pattern << " offset = " << patternOffset
0167                                 << " phi_local = " << data.phi_local;
0168
0169   /// Local Phi Bend is always zero. Until we start using it.
0170   data.phi_bend_local = 0;
0171
0172   return data;  //return LUT result
0173 }
0174
0175 void CSCSectorReceiverLUT::fillLocalPhiLUT() {
0176   // read data in from a file... Add this later.
0177 }
0178
0179 lclphidat CSCSectorReceiverLUT::localPhi(
0180     const int strip, const int pattern, const int quality, const int lr, const bool gangedME1a) const {
0181   lclphiadd theadd;
0182
0183   theadd.strip = strip;
0184   theadd.clct_pattern = pattern & 0x7;
0185   theadd.pattern_type = (pattern & 0x8) >> 3;
0186   theadd.quality = quality;
0187   theadd.lr = lr;
0188   theadd.spare = 0;
0189
0190   return localPhi(theadd, gangedME1a);
0191 }
0192
0193 lclphidat CSCSectorReceiverLUT::localPhi(unsigned address, const bool gangedME1a) const {
0194   lclphidat result;
0195   lclphiadd theadd(address);
0196
0197   if (useMiniLUTs && isTMB07) {
0198     result = CSCSectorReceiverMiniLUT::calcLocalPhiMini(address, gangedME1a);
0199   } else if (LUTsFromFile)
0200     result = me_lcl_phi[address];
0201   else
0202     result = calcLocalPhi(theadd);
0203
0204   return result;
0205 }
0206
0207 lclphidat CSCSectorReceiverLUT::localPhi(lclphiadd address, const bool gangedME1a) const {
0208   lclphidat result;
0209
0210   if (useMiniLUTs && isTMB07) {
0211     result = CSCSectorReceiverMiniLUT::calcLocalPhiMini(address.toint(), gangedME1a);
0212   } else if (LUTsFromFile)
0213     result = me_lcl_phi[address.toint()];
0214   else
0215     result = calcLocalPhi(address);
0216
0217   return result;
0218 }
0219
0220 double CSCSectorReceiverLUT::getGlobalPhiValue(const CSCLayer* thelayer,
0221                                                const unsigned& strip,
0222                                                const unsigned& wire_group) const {
0223   double result = 0.0;
0224   //CSCLayerGeometry* thegeom;
0225   //LocalPoint lp;
0226   //GlobalPoint gp;
0227
0228   try {
0229     //thegeom = const_cast<CSCLayerGeometry*>(thelayer->geometry());
0230     //lp = thegeom->stripWireGroupIntersection(strip, wire_group);
0231     //gp = thelayer->surface().toGlobal(lp);
0232     result = thelayer->centerOfStrip(strip).phi();  //gp.phi();
0233
0234     if (result < 0.)
0235       result += 2. * M_PI;
0236   } catch (edm::Exception& e) {
0237     LogDebug("CSCSectorReceiverLUT|getGlobalPhiValue") << e.what();
0238   }
0239
0240   return result;
0241 }
0242
0243 gblphidat CSCSectorReceiverLUT::calcGlobalPhiME(const gblphiadd& address) const {
0244   gblphidat result(0);
0245   const CSCChamber* thechamber = nullptr;
0246   const CSCLayer* thelayer = nullptr;
0247   const CSCLayerGeometry* layergeom = nullptr;
0248   int cscid = address.cscid;
0249   unsigned wire_group = address.wire_group;
0250   unsigned local_phi = address.phi_local;
0251   const double sectorOffset =
0252       (CSCTFConstants::SECTOR1_CENT_RAD - CSCTFConstants::SECTOR_RAD / 2.) + (_sector - 1) * M_PI / 3.;
0253
0254   //Number of global phi units per radian.
0255   constexpr int maxPhiG = 1 << CSCBitWidths::kGlobalPhiDataBitWidth;
0256   double binPhiG = static_cast<double>(maxPhiG) / CSCTFConstants::SECTOR_RAD;
0257
0258   // We will use these to convert the local phi into radians.
0259   constexpr unsigned int maxPhiL = 1 << CSCBitWidths::kLocalPhiDataBitWidth;
0260   const double binPhiL = static_cast<double>(maxPhiL) / (2. * CSCConstants::MAX_NUM_STRIPS_RUN1);
0261
0262   if (cscid < CSCTriggerNumbering::minTriggerCscId()) {
0263     edm::LogWarning("CSCSectorReceiverLUT|getGlobalPhiValue")
0264         << " warning: cscId " << cscid << " is out of bounds [" << CSCTriggerNumbering::minTriggerCscId() << "-"
0265         << CSCTriggerNumbering::maxTriggerCscId() << "]\n";
0266     throw cms::Exception("CSCSectorReceiverLUT")
0267         << "+++ Value of CSC ID, " << cscid << ", is out of bounds [" << CSCTriggerNumbering::minTriggerCscId() << "-"
0268         << CSCTriggerNumbering::maxTriggerCscId() << "] +++\n";
0269   }
0270
0271   if (cscid > CSCTriggerNumbering::maxTriggerCscId()) {
0272     edm::LogWarning("CSCSectorReceiverLUT|getGlobalPhiValue")
0273         << " warning: cscId " << cscid << " is out of bounds [" << CSCTriggerNumbering::minTriggerCscId() << "-"
0274         << CSCTriggerNumbering::maxTriggerCscId() << "]\n";
0275     throw cms::Exception("CSCSectorReceiverLUT")
0276         << "+++ Value of CSC ID, " << cscid << ", is out of bounds [" << CSCTriggerNumbering::minTriggerCscId() << "-"
0277         << CSCTriggerNumbering::maxTriggerCscId() << "] +++\n";
0278   }
0279
0280   if (wire_group >= 1 << 5) {
0281     edm::LogWarning("CSCSectorReceiverLUT|getGlobalPhiValue")
0282         << "warning: wire_group" << wire_group << " is out of bounds (1-" << ((1 << 5) - 1) << "]\n";
0283     throw cms::Exception("CSCSectorReceiverLUT")
0284         << "+++ Value of wire_group, " << wire_group << ", is out of bounds (1-" << ((1 << 5) - 1) << "] +++\n";
0285   }
0286
0287   if (local_phi >= maxPhiL) {
0288     edm::LogWarning("CSCSectorReceiverLUT|getGlobalPhiValue")
0289         << "warning: local_phi" << local_phi << " is out of bounds [0-" << maxPhiL << ")\n";
0290     throw cms::Exception("CSCSectorReceiverLUT")
0291         << "+++ Value of local_phi, " << local_phi << ", is out of bounds [0-, " << maxPhiL << ") +++\n";
0292   }
0293
0294   try {
0295     int ring = CSCTriggerNumbering::ringFromTriggerLabels(_station, cscid);
0296     int chid = CSCTriggerNumbering::chamberFromTriggerLabels(_sector, _subsector, _station, cscid);
0297     CSCDetId detid(_endcap, _station, ring, chid, 0);
0298     thechamber = const_cast<const CSCChamber*>(csc_g->chamber(detid));
0299     if (thechamber) {
0300       layergeom = thechamber->layer(CSCConstants::KEY_CLCT_LAYER)->geometry();
0301       thelayer = thechamber->layer(CSCConstants::KEY_CLCT_LAYER);
0302       const int nStrips = layergeom->numberOfStrips();
0303       // PhiL is the strip number converted into some units between 0 and
0304       // 1023.  When we did the conversion in fillLocalPhiTable(), we did
0305       // not know for which chamber we do it (and, therefore, how many strips
0306       // it has), and always used the maximum possible number of strips
0307       // per chamber, MAX_NUM_STRIPS_RUN1=80.  Now, since we know the chamber id
0308       // and how many strips the chamber has, we can re-adjust the scale.
0309       //const double scale = static_cast<double>(CSCConstants::MAX_NUM_STRIPS_RUN1)/nStrips;
0310
0311       int strip = 0, halfstrip = 0;
0312
0313       halfstrip = static_cast<int>(local_phi / binPhiL);
0314       strip = halfstrip / 2;
0315
0316       // Find the phi width of the chamber and the position of its "left"
0317       // (lower phi) edge (both in radians).
0318       // Phi positions of the centers of the first and of the last strips
0319       // in the chamber.
0320       const double phi_f = getGlobalPhiValue(thelayer, 1, wire_group);
0321       const double phi_l = getGlobalPhiValue(thelayer, nStrips, wire_group);
0322       // Phi widths of the half-strips at both ends of the chamber;
0323       // surprisingly, they are not the same.
0324       const double hsWidth_f = fabs(getGlobalPhiValue(thelayer, 2, wire_group) - phi_f) / 2.;
0325       const double hsWidth_l = fabs(phi_l - getGlobalPhiValue(thelayer, nStrips - 1, wire_group)) / 2.;
0326
0327       // The "natural" match between the strips and phi values -- when
0328       // a larger strip number corresponds to a larger phi value, i.e. strips
0329       // are counted clockwise if we look at them from the inside of the
0330       // detector -- is reversed for some stations.  At the moment, these
0331       // are stations 3 and 4 of the 1st endcap, and stations 1 and 2 of
0332       // the 2nd endcap.  Instead of using
0333       // if ((theEndcap == 1 && theStation <= 2) ||
0334       // (theEndcap == 2 && theStation >= 3)),
0335       // we get the order from the phi values of the first and the last strip
0336       // in a chamber, just in case the counting scheme changes in the future.
0337       // Once we know how the strips are counted, we can go from the middle
0338       // of the strips to their outer edges.
0339       bool clockwiseOrder;
0340       double leftEdge, rightEdge;
0341       if (fabs(phi_f - phi_l) < M_PI) {
0342         if (phi_f < phi_l)
0343           clockwiseOrder = true;
0344         else
0345           clockwiseOrder = false;
0346       } else {  // the chamber crosses the phi = pi boundary
0347         if (phi_f < phi_l)
0348           clockwiseOrder = false;
0349         else
0350           clockwiseOrder = true;
0351       }
0352       if (clockwiseOrder) {
0353         leftEdge = phi_f - hsWidth_f;
0354         rightEdge = phi_l + hsWidth_l;
0355       } else {
0356         leftEdge = phi_l - hsWidth_l;
0357         rightEdge = phi_f + hsWidth_f;
0358       }
0359       if (fabs(phi_f - phi_l) >= M_PI) {
0360         rightEdge += 2. * M_PI;
0361       }
0362       //double chamberWidth = (rightEdge - leftEdge);
0363
0364       // Chamber offset, relative to the edge of the sector.
0365       //double chamberOffset = leftEdge - sectorOffset;
0366       //if (chamberOffset < -M_PI) chamberOffset += 2*M_PI;
0367
0368       double temp_phi = 0.0, strip_phi = 0.0, delta_phi = 0.0;
0369       double distFromHalfStripCenter = 0.0, halfstripWidth = 0.0;
0370
0371       if (strip < nStrips) {
0372         // Approximate distance from the center of the half-strip to the center
0373         // of this phil bin, in units of half-strip width.
0374         distFromHalfStripCenter = (local_phi + 0.5) / binPhiL - halfstrip - 0.5;
0375         // Half-strip width (in rad), calculated as the half-distance between
0376         // the adjacent strips.  Since in the current ORCA implementation
0377         // the half-strip width changes from strip to strip, base the choice
0378         // of the adjacent strip on the half-strip number.
0379         if ((halfstrip % 2 == 0 && halfstrip != 0) || halfstrip == 2 * nStrips - 1) {
0380           halfstripWidth = fabs(getGlobalPhiValue(thelayer, strip + 1, wire_group) -
0381                                 getGlobalPhiValue(thelayer, strip, wire_group)) /
0382                            2.;
0383         } else {
0384           halfstripWidth = fabs(getGlobalPhiValue(thelayer, strip + 1, wire_group) -
0385                                 getGlobalPhiValue(thelayer, strip + 2, wire_group)) /
0386                            2.;
0387         }
0388         // Correction for the strips crossing the 180 degree boundary.
0389         if (halfstripWidth > M_PI / 2.)
0390           halfstripWidth = M_PI - halfstripWidth;
0391         // Phi at the center of the strip.
0392         strip_phi = getGlobalPhiValue(thelayer, strip + 1, wire_group);
0393         // Distance between the center of the strip and the phil position.
0394         delta_phi = halfstripWidth * (((halfstrip % 2) - 0.5) + distFromHalfStripCenter);
0395         if (clockwiseOrder)
0396           temp_phi = strip_phi + delta_phi;
0397         else
0398           temp_phi = strip_phi - delta_phi;
0399       } else {
0400         // PhiL values that do not have corresponding strips (the chamber
0401         // has less than 80 strips assumed in fillLocalPhi).  It does not
0402         // really matter what we do with these values; at the moment, just
0403         // set them to the phis of the edges of the chamber.
0404         if (clockwiseOrder)
0405           temp_phi = rightEdge;
0406         else
0407           temp_phi = leftEdge;
0408       }
0409
0410       // Finally, subtract the sector offset and convert to the scale of
0411       // the global phi.
0412
0413       temp_phi -= sectorOffset;
0414
0415       if (temp_phi < 0.)
0416         temp_phi += 2. * M_PI;
0417
0418       temp_phi *= binPhiG;
0419
0420       if (temp_phi < 0.) {
0421         result.global_phi = 0;
0422       } else if (temp_phi >= maxPhiG) {
0423         result.global_phi = maxPhiG - 1;
0424       } else {
0425         result.global_phi = static_cast<unsigned short>(temp_phi);
0426       }
0427
0428       LogDebug("CSCSectorReceiverLUT") << "local_phi = " << local_phi << " halfstrip = " << halfstrip
0429                                        << " strip = " << strip
0430                                        << " distFromHalfStripCenter = " << distFromHalfStripCenter
0431                                        << " halfstripWidth = " << halfstripWidth
0432                                        << " strip phi = " << strip_phi / (M_PI / 180.)
0433                                        << " temp_phi = " << temp_phi * CSCTFConstants::SECTOR_DEG / maxPhiG
0434                                        << " global_phi = " << result.global_phi << " "
0435                                        << result.global_phi * CSCTFConstants::SECTOR_DEG / maxPhiG;
0436     }
0437   } catch (edm::Exception& e) {
0438     edm::LogError("CSCSectorReceiverLUT|getGlobalPhiValue") << e.what();
0439   }
0440
0441   return result;
0442 }
0443
0444 gblphidat CSCSectorReceiverLUT::globalPhiME(int phi_local, int wire_group, int cscid, const bool gangedME1a) const {
0445   gblphidat result;
0446   gblphiadd theadd;
0447   theadd.phi_local = phi_local;
0448   theadd.wire_group = ((1 << 5) - 1) & (wire_group >> 2);  // want 2-7 of wg
0449   theadd.cscid = cscid;
0450
0451   if (useMiniLUTs && isTMB07)
0452     result = CSCSectorReceiverMiniLUT::calcGlobalPhiMEMini(
0453         _endcap, _sector, _station, _subsector, theadd.toint(), gangedME1a);
0454   else if (LUTsFromFile)
0455     result = me_global_phi[theadd.toint()];
0456   else
0457     result = calcGlobalPhiME(theadd);
0458
0459   return result;
0460 }
0461
0462 gblphidat CSCSectorReceiverLUT::globalPhiME(unsigned address, const bool gangedME1a) const {
0463   gblphidat result;
0464
0465   if (useMiniLUTs && isTMB07)
0466     result = CSCSectorReceiverMiniLUT::calcGlobalPhiMEMini(_endcap, _sector, _station, _subsector, address, gangedME1a);
0467   else if (LUTsFromFile)
0468     result = me_global_phi[address];
0469   else
0470     result = calcGlobalPhiME(gblphiadd(address));
0471
0472   return result;
0473 }
0474
0475 gblphidat CSCSectorReceiverLUT::globalPhiME(gblphiadd address, const bool gangedME1a) const {
0476   gblphidat result;
0477
0478   if (useMiniLUTs && isTMB07)
0479     result = CSCSectorReceiverMiniLUT::calcGlobalPhiMEMini(
0480         _endcap, _sector, _station, _subsector, address.toint(), gangedME1a);
0481   else if (LUTsFromFile)
0482     result = me_global_phi[address.toint()];
0483   else
0484     result = calcGlobalPhiME(address);
0485
0486   return result;
0487 }
0488
0489 gblphidat CSCSectorReceiverLUT::calcGlobalPhiMB(const gblphidat& csclut) const {
0490   gblphidat dtlut;
0491
0492   // The following method was ripped from D. Holmes' LUT conversion program
0493   // modifications from Darin and GP
0494   int GlobalPhiMin = (_subsector == 1) ? 0x42 : 0x800;   // (0.999023 : 31 in degrees)
0495   int GlobalPhiMax = (_subsector == 1) ? 0x7ff : 0xfbd;  // (30.985 : 60.986 in degrees)
0496   double GlobalPhiShift = (1.0 * GlobalPhiMin + (GlobalPhiMax - GlobalPhiMin) / 2.0);
0497
0498   double dt_out = static_cast<double>(csclut.global_phi) - GlobalPhiShift;
0499
0500   // these numbers are 62 deg / 1 rad (CSC phi scale vs. DT phi scale)
0501   dt_out = (dt_out / 1982) * 2145;  //CSC phi 62 degrees; DT phi 57.3 degrees
0502
0503   if (dt_out >= 0)  // msb != 1
0504   {
0505     dtlut.global_phi = 0x7ff & static_cast<unsigned>(dt_out);
0506   } else {
0507     dtlut.global_phi = static_cast<unsigned>(-dt_out);
0508     dtlut.global_phi = ~dtlut.global_phi;
0509     dtlut.global_phi |= 0x800;
0510   }
0511
0512   return dtlut;
0513 }
0514
0515 gblphidat CSCSectorReceiverLUT::globalPhiMB(int phi_local, int wire_group, int cscid, const bool gangedME1a) const {
0516   gblphiadd address;
0517   gblphidat result;
0518
0519   address.cscid = cscid;
0520   address.wire_group = ((1 << 5) - 1) & (wire_group >> 2);
0521   address.phi_local = phi_local;
0522
0523   // comment for now
0524   //  if(useMiniLUTs && isTMB07) result = CSCSectorReceiverMiniLUT::calcGlobalPhiMBMini(_endcap, _sector, _subsector, address.toint());
0525   //else
0526   if (LUTsFromFile)
0527     result = mb_global_phi[address.toint()];
0528   else
0529     result = calcGlobalPhiMB(globalPhiME(address, gangedME1a));
0530
0531   return result;
0532 }
0533
0534 gblphidat CSCSectorReceiverLUT::globalPhiMB(unsigned address, const bool gangedME1a) const {
0535   gblphidat result;
0536   gblphiadd theadd(address);
0537
0538   //if(useMiniLUTs && isTMB07) result = CSCSectorReceiverMiniLUT::calcGlobalPhiMBMini(_endcap, _sector, _subsector, address);
0539   //else
0540   if (LUTsFromFile)
0541     result = mb_global_phi[theadd.toint()];
0542   else
0543     result = calcGlobalPhiMB(globalPhiME(address, gangedME1a));
0544
0545   return result;
0546 }
0547
0548 gblphidat CSCSectorReceiverLUT::globalPhiMB(gblphiadd address, const bool gangedME1a) const {
0549   gblphidat result;
0550
0551   //if(useMiniLUTs && isTMB07) result = CSCSectorReceiverMiniLUT::calcGlobalPhiMBMini(_endcap, _sector, _subsector, address.toint());
0552   //else
0553   if (LUTsFromFile)
0554     result = mb_global_phi[address.toint()];
0555   else
0556     result = calcGlobalPhiMB(globalPhiME(address, gangedME1a));
0557
0558   return result;
0559 }
0560
0561 double CSCSectorReceiverLUT::getGlobalEtaValue(const unsigned& thecscid,
0562                                                const unsigned& thewire_group,
0563                                                const unsigned& thephi_local) const {
0564   double result = 0.0;
0565   unsigned wire_group = thewire_group;
0566   int cscid = thecscid;
0567   unsigned phi_local = thephi_local;
0568
0569   // Flag to be set if one wants to apply phi corrections ONLY in ME1/1.
0570   // Turn it into a parameter?
0571   bool me1ir_only = false;
0572
0573   if (cscid < CSCTriggerNumbering::minTriggerCscId() || cscid > CSCTriggerNumbering::maxTriggerCscId()) {
0574     edm::LogWarning("CSCSectorReceiverLUT|getEtaValue")
0575         << " warning: cscId " << cscid << " is out of bounds [1-" << CSCTriggerNumbering::maxTriggerCscId() << "]\n";
0576     cscid = CSCTriggerNumbering::maxTriggerCscId();
0577   }
0578
0579   CSCLayerGeometry* layerGeom = nullptr;
0580   const unsigned numBins = 1 << 2;  // 4 local phi bins
0581
0582   if (phi_local > numBins - 1) {
0583     edm::LogWarning("CSCSectorReceiverLUT|getEtaValue")
0584         << "warning: phiL " << phi_local << " is out of bounds [0-" << numBins - 1 << "]\n";
0585     phi_local = numBins - 1;
0586   }
0587   try {
0588     int ring = CSCTriggerNumbering::ringFromTriggerLabels(_station, cscid);
0589     int chid = CSCTriggerNumbering::chamberFromTriggerLabels(_sector, _subsector, _station, cscid);
0590     CSCDetId detid(_endcap, _station, ring, chid, 0);
0591     const CSCChamber* thechamber = const_cast<const CSCChamber*>(csc_g->chamber(detid));
0592     if (thechamber) {
0593       layerGeom = const_cast<CSCLayerGeometry*>(thechamber->layer(CSCConstants::KEY_ALCT_LAYER)->geometry());
0594       const unsigned nWireGroups = layerGeom->numberOfWireGroups();
0595
0596       // Check wire group numbers; expect them to be counted from 0, as in
0597       // CorrelatedLCTDigi class.
0598       if (wire_group >= nWireGroups) {
0599         edm::LogWarning("CSCSectorReceiverLUT|getEtaValue")
0600             << "warning: wireGroup " << wire_group << " is out of bounds [0-" << nWireGroups << ")\n";
0601         wire_group = nWireGroups - 1;
0602       }
0603       // Convert to [1; nWireGroups] range used in geometry methods.
0604       wire_group += 1;
0605
0606       // If me1ir_only is set, apply phi corrections only in ME1/1.
0607       if (me1ir_only && (_station != 1 || CSCTriggerNumbering::ringFromTriggerLabels(_station, cscid) != 1)) {
0608         result = thechamber->layer(CSCConstants::KEY_ALCT_LAYER)->centerOfWireGroup(wire_group).eta();
0609       } else {
0610         const unsigned nStrips = layerGeom->numberOfStrips();
0611         const unsigned nStripsPerBin = CSCConstants::MAX_NUM_STRIPS_RUN1 / numBins;
0612         /**
0613        * Calculate Eta correction
0614        */
0615
0616         // Check that no strips will be left out.
0617         if (nStrips % numBins != 0 || CSCConstants::MAX_NUM_STRIPS_RUN1 % numBins != 0)
0618           edm::LogWarning("CSCSectorReceiverLUT")
0619               << "getGlobalEtaValue warning: number of strips " << nStrips << " (" << CSCConstants::MAX_NUM_STRIPS_RUN1
0620               << ") is not divisible by numBins " << numBins << " Station " << _station << " sector " << _sector
0621               << " subsector " << _subsector << " cscid " << cscid << "\n";
0622
0623         unsigned maxStripPrevBin = 0, maxStripThisBin = 0;
0624         unsigned correctionStrip;
0625         LocalPoint lPoint;
0626         GlobalPoint gPoint;
0627         // Bins phi_local and find the the middle strip for each bin.
0628         maxStripThisBin = nStripsPerBin * (phi_local + 1);
0629         if (maxStripThisBin <= nStrips) {
0630           correctionStrip = nStripsPerBin / 2 * (2 * phi_local + 1);
0631         } else {
0632           // If the actual number of strips in the chamber is smaller than
0633           // the number of strips corresponding to the right edge of this phi
0634           // local bin, we take the middle strip between number of strips
0635           // at the left edge of the bin and the actual number of strips.
0636           maxStripPrevBin = nStripsPerBin * phi_local;
0637           correctionStrip = (nStrips + maxStripPrevBin) / 2;
0638         }
0639
0640         lPoint = layerGeom->stripWireGroupIntersection(correctionStrip, wire_group);
0641         gPoint = thechamber->layer(CSCConstants::KEY_ALCT_LAYER)->surface().toGlobal(lPoint);
0642
0643         // end calc of eta correction.
0644         result = gPoint.eta();
0645       }
0646     }
0647   } catch (cms::Exception& e) {
0648     LogDebug("CSCSectorReceiver|OutofBoundInput") << e.what();
0649   }
0650
0651   return std::fabs(result);
0652 }
0653
0654 gbletadat CSCSectorReceiverLUT::calcGlobalEtaME(const gbletaadd& address) const {
0655   gbletadat result;
0656   double float_eta = getGlobalEtaValue(address.cscid, address.wire_group, address.phi_local);
0657   unsigned int_eta = 0;
0658   unsigned bend_global = 0;  // not filled yet... will change when it is.
0659   const double etaPerBin = (CSCTFConstants::maxEta - CSCTFConstants::minEta) / CSCTFConstants::etaBins;
0660   const unsigned me12EtaCut = 56;
0661
0662   if ((float_eta < CSCTFConstants::minEta) || (float_eta >= CSCTFConstants::maxEta)) {
0663     edm::LogWarning("CSCSectorReceiverLUT:OutOfBounds")
0664         << "CSCSectorReceiverLUT warning: float_eta = " << float_eta << " minEta = " << CSCTFConstants::minEta
0665         << " maxEta = " << CSCTFConstants::maxEta << "   station " << _station << " sector " << _sector << " chamber "
0666         << address.cscid << " wire group " << address.wire_group;
0667
0668     throw cms::Exception("CSCSectorReceiverLUT")
0669         << "+++ Value of CSC ID, " << float_eta << ", is out of bounds [" << CSCTFConstants::minEta << "-"
0670         << CSCTFConstants::maxEta << ") +++\n";
0671
0672     //if (float_eta < CSCTFConstants::minEta)
0673     //result.global_eta = 0;
0674     //else if (float_eta >= CSCTFConstants::maxEta)
0675     //result.global_eta = CSCTFConstants::etaBins - 1;
0676   } else {
0677     float_eta -= CSCTFConstants::minEta;
0678     float_eta = float_eta / etaPerBin;
0679     int_eta = static_cast<unsigned>(float_eta);
0680     /* Commented until I find out its use.
0681       // Fine-tune eta boundary between DT and CSC.
0682       if ((intEta == L1MuCSCSetup::CscEtaStart() && (L1MuCSCSetup::CscEtaStartCorr() > 0.) ) ||
0683       (intEta == L1MuCSCSetup::CscEtaStart() - 1 && (L1MuCSCSetup::CscEtaStartCorr() < 0.) ) ) {
0684     bitEta = (thisEta-minEta-L1MuCSCSetup::CscEtaStartCorr())/EtaPerBin;
0685     intEta = static_cast<int>(bitEta);
0686       }
0687       */
0688     if (_station == 1 && address.cscid >= static_cast<unsigned>(CSCTriggerNumbering::minTriggerCscId()) &&
0689         address.cscid <= static_cast<unsigned>(CSCTriggerNumbering::maxTriggerCscId())) {
0690       unsigned ring = CSCTriggerNumbering::ringFromTriggerLabels(_station, address.cscid);
0691
0692       if (ring == 1 && int_eta < me12EtaCut) {
0693         int_eta = me12EtaCut;
0694       } else if (ring == 2 && int_eta >= me12EtaCut) {
0695         int_eta = me12EtaCut - 1;
0696       }
0697     }
0698     result.global_eta = int_eta;
0699   }
0700   result.global_bend = bend_global;
0701
0702   return result;
0703 }
0704
0705 gbletadat CSCSectorReceiverLUT::globalEtaME(
0706     int tphi_bend, int tphi_local, int twire_group, int tcscid, const bool gangedME1a) const {
0707   gbletadat result;
0708   gbletaadd theadd;
0709
0710   theadd.phi_bend = tphi_bend;
0711   theadd.phi_local = (tphi_local >> (CSCBitWidths::kLocalPhiDataBitWidth - 2)) & 0x3;  // want 2 msb of local phi
0712   theadd.wire_group = twire_group;
0713   theadd.cscid = tcscid;
0714   if (useMiniLUTs && isTMB07)
0715     result = CSCSectorReceiverMiniLUT::calcGlobalEtaMEMini(
0716         _endcap, _sector, _station, _subsector, theadd.toint(), gangedME1a);
0717   else if (LUTsFromFile)
0718     result = me_global_eta[theadd.toint()];
0719   else
0720     result = calcGlobalEtaME(theadd);
0721
0722   return result;
0723 }
0724
0725 gbletadat CSCSectorReceiverLUT::globalEtaME(unsigned address, const bool gangedME1a) const {
0726   gbletadat result;
0727   gbletaadd theadd(address);
0728
0729   if (useMiniLUTs && isTMB07)
0730     result = CSCSectorReceiverMiniLUT::calcGlobalEtaMEMini(_endcap, _sector, _station, _subsector, address, gangedME1a);
0731   else if (LUTsFromFile)
0732     result = me_global_eta[address];
0733   else
0734     result = calcGlobalEtaME(theadd);
0735   return result;
0736 }
0737
0738 gbletadat CSCSectorReceiverLUT::globalEtaME(gbletaadd address, const bool gangedME1a) const {
0739   gbletadat result;
0740
0741   if (useMiniLUTs && isTMB07)
0742     result = CSCSectorReceiverMiniLUT::calcGlobalEtaMEMini(
0743         _endcap, _sector, _station, _subsector, address.toint(), gangedME1a);
0744   else if (LUTsFromFile)
0745     result = me_global_eta[address.toint()];
0746   else
0747     result = calcGlobalEtaME(address);
0748   return result;
0749 }
0750
0751 std::string CSCSectorReceiverLUT::encodeFileIndex() const {
0752   std::string fileName = "";
0753   if (_station == 1) {
0754     if (_subsector == 1)
0755       fileName += "1a";
0756     if (_subsector == 2)
0757       fileName += "1b";
0758   } else if (_station == 2)
0759     fileName += "2";
0760   else if (_station == 3)
0761     fileName += "3";
0762   else if (_station == 4)
0763     fileName += "4";
0764   fileName += "End";
0765   if (_endcap == 1)
0766     fileName += "1";
0767   else
0768     fileName += "2";
0769   fileName += "Sec";
0770   if (_sector == 1)
0771     fileName += "1";
0772   else if (_sector == 2)
0773     fileName += "2";
0774   else if (_sector == 3)
0775     fileName += "3";
0776   else if (_sector == 4)
0777     fileName += "4";
0778   else if (_sector == 5)
0779     fileName += "5";
0780   else if (_sector == 6)
0781     fileName += "6";
0782   fileName += "LUT";
0783   return fileName;
0784 }
0785
0786 void CSCSectorReceiverLUT::readLUTsFromFile() {
0787   if (!me_lcl_phi_loaded) {
0788     me_lcl_phi = new lclphidat[1 << CSCBitWidths::kLocalPhiAddressWidth];
0789     std::string fName(me_lcl_phi_file.fullPath());
0790     std::ifstream LocalPhiLUT;
0791
0792     edm::LogInfo("CSCSectorReceiverLUT") << "Loading SR LUT: " << fName;
0793
0794     if (isBinary) {
0795       LocalPhiLUT.open(fName.c_str(), std::ios::binary);
0796       LocalPhiLUT.seekg(0, std::ios::end);
0797       int length = LocalPhiLUT.tellg();
0798       if (length == (1 << CSCBitWidths::kLocalPhiAddressWidth) * sizeof(short)) {
0799         LocalPhiLUT.seekg(0, std::ios::beg);
0800         LocalPhiLUT.read(reinterpret_cast<char*>(me_lcl_phi), length);
0801         LocalPhiLUT.close();
0802       } else
0803         edm::LogError("CSCSectorReceiverLUT") << "File " << fName << " is incorrect size!";
0804       LocalPhiLUT.close();
0805     } else {
0806       LocalPhiLUT.open(fName.c_str());
0807       unsigned i = 0;
0808       unsigned short temp = 0;
0809       while (!LocalPhiLUT.eof() && i < 1 << CSCBitWidths::kLocalPhiAddressWidth) {
0810         LocalPhiLUT >> temp;
0811         me_lcl_phi[i++] = temp;
0812       }
0813       LocalPhiLUT.close();
0814     }
0815   }
0816   if (!me_global_phi) {
0817     me_global_phi = new gblphidat[1 << CSCBitWidths::kGlobalPhiAddressWidth];
0818     std::string fName = me_gbl_phi_file.fullPath();
0819     std::ifstream GlobalPhiLUT;
0820
0821     edm::LogInfo("CSCSectorReceiverLUT") << "Loading SR LUT: " << fName;
0822
0823     if (isBinary) {
0824       GlobalPhiLUT.open(fName.c_str(), std::ios::binary);
0825       GlobalPhiLUT.seekg(0, std::ios::end);
0826       int length = GlobalPhiLUT.tellg();
0827       if (length == (1 << CSCBitWidths::kGlobalPhiAddressWidth) * sizeof(short)) {
0828         GlobalPhiLUT.seekg(0, std::ios::beg);
0829         GlobalPhiLUT.read(reinterpret_cast<char*>(me_global_phi), length);
0830       } else
0831         edm::LogError("CSCSectorReceiverLUT") << "File " << fName << " is incorrect size!";
0832       GlobalPhiLUT.close();
0833     } else {
0834       GlobalPhiLUT.open(fName.c_str());
0835       unsigned short temp = 0;
0836       unsigned i = 0;
0837       while (!GlobalPhiLUT.eof() && i < 1 << CSCBitWidths::kGlobalPhiAddressWidth) {
0838         GlobalPhiLUT >> temp;
0839         me_global_phi[i++] = temp;
0840       }
0841       GlobalPhiLUT.close();
0842     }
0843   }
0844   if (!mb_global_phi && _station == 1)  // MB lut only in station one.
0845   {
0846     mb_global_phi = new gblphidat[1 << CSCBitWidths::kGlobalPhiAddressWidth];
0847     std::string fName = mb_gbl_phi_file.fullPath();
0848     std::ifstream GlobalPhiLUT;
0849
0850     edm::LogInfo("CSCSectorReceiverLUT") << "Loading SR LUT: " << fName;
0851
0852     if (isBinary) {
0853       GlobalPhiLUT.open(fName.c_str(), std::ios::binary);
0854       GlobalPhiLUT.seekg(0, std::ios::end);
0855       int length = GlobalPhiLUT.tellg();
0856       if (length == (1 << CSCBitWidths::kGlobalPhiAddressWidth) * sizeof(short)) {
0857         GlobalPhiLUT.seekg(0, std::ios::beg);
0858         GlobalPhiLUT.read(reinterpret_cast<char*>(mb_global_phi), length);
0859       } else
0860         edm::LogError("CSCSectorReceiverLUT") << "File " << fName << " is incorrect size!";
0861       GlobalPhiLUT.close();
0862     } else {
0863       GlobalPhiLUT.open(fName.c_str());
0864       unsigned short temp = 0;
0865       unsigned i = 0;
0866       while (!GlobalPhiLUT.eof() && i < 1 << CSCBitWidths::kGlobalPhiAddressWidth) {
0867         GlobalPhiLUT >> temp;
0868         mb_global_phi[i++] = temp;
0869       }
0870       GlobalPhiLUT.close();
0871     }
0872   }
0873   if (!me_global_eta) {
0874     me_global_eta = new gbletadat[1 << CSCBitWidths::kGlobalEtaAddressWidth];
0875     std::string fName = me_gbl_eta_file.fullPath();
0876     std::ifstream GlobalEtaLUT;
0877
0878     edm::LogInfo("CSCSectorReceiverLUT") << "Loading SR LUT: " << fName;
0879
0880     if (isBinary) {
0881       GlobalEtaLUT.open(fName.c_str(), std::ios::binary);
0882       GlobalEtaLUT.seekg(0, std::ios::end);
0883       int length = GlobalEtaLUT.tellg();
0884       if (length == (1 << CSCBitWidths::kGlobalEtaAddressWidth) * sizeof(short)) {
0885         GlobalEtaLUT.seekg(0, std::ios::beg);
0886         GlobalEtaLUT.read(reinterpret_cast<char*>(me_global_eta), length);
0887       } else
0888         edm::LogError("CSCSectorReceiverLUT") << "File " << fName << " is incorrect size!";
0889       GlobalEtaLUT.close();
0890     } else {
0891       GlobalEtaLUT.open(fName.c_str());
0892       unsigned short temp = 0;
0893       unsigned i = 0;
0894       while (!GlobalEtaLUT.eof() && i < 1 << CSCBitWidths::kGlobalEtaAddressWidth) {
0895         GlobalEtaLUT >> temp;
0896         me_global_eta[i++] = temp;
0897       }
0898       GlobalEtaLUT.close();
0899     }
0900   }
0901 }