Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-04-21 02:53:23

0001 #include <cmath>
0002 #include <memory>
0003 #include <vector>
0004 #include <fstream>
0005 #include <iostream>
0006 
0007 #include "TFile.h"
0008 #include "TTree.h"
0009 
0010 #include "FWCore/Framework/interface/Frameworkfwd.h"
0011 #include "FWCore/Framework/interface/EDAnalyzer.h"
0012 #include "FWCore/Framework/interface/Event.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0015 #include "FWCore/Framework/interface/EventSetup.h"
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017 
0018 #include "DataFormats/GeometryVector/interface/LocalPoint.h"
0019 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0020 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0021 #include "DataFormats/MuonDetId/interface/CSCTriggerNumbering.h"
0022 #include "DataFormats/CSCDigi/interface/CSCConstants.h"
0023 
0024 #include "Geometry/CSCGeometry/interface/CSCGeometry.h"
0025 #include "Geometry/CSCGeometry/interface/CSCChamber.h"
0026 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0027 #include "Geometry/CSCGeometry/interface/CSCLayerGeometry.h"
0028 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0029 
0030 #include "helper.h"
0031 
0032 class MakeCoordLUT : public edm::EDAnalyzer {
0033 public:
0034   explicit MakeCoordLUT(const edm::ParameterSet&);
0035   virtual ~MakeCoordLUT();
0036 
0037 private:
0038   //virtual void beginJob();
0039   //virtual void endJob();
0040 
0041   virtual void beginRun(const edm::Run&, const edm::EventSetup&);
0042   virtual void endRun(const edm::Run&, const edm::EventSetup&);
0043 
0044   virtual void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup);
0045 
0046   // Generate LUTs
0047   void generateLUTs();
0048   void generateLUTs_init();
0049   void generateLUTs_run();
0050   void generateLUTs_final();
0051 
0052   // Validate LUTs
0053   void validateLUTs();
0054 
0055   // Write LUT files
0056   void writeFiles();
0057 
0058   // Construct CSCDetId
0059   CSCDetId getCSCDetId(int endcap, int sector, int subsector, int station, int cscid, bool isME1A) const;
0060 
0061   // Is strip phi counter-clockwise
0062   bool isStripPhiCounterClockwise(const CSCDetId& cscDetId) const;
0063 
0064   // Get strip pitch
0065   double getStripPitch(const CSCDetId& cscDetId) const;
0066 
0067   // Get global phi in degrees
0068   double getGlobalPhi(
0069       int endcap, int sector, int subsector, int station, int cscid, bool isME1A, int wiregroup, int halfstrip) const;
0070   double getGlobalPhiFullstrip(
0071       int endcap, int sector, int subsector, int station, int cscid, bool isME1A, int wiregroup, int fullstrip) const;
0072 
0073   // Get global theta in degrees
0074   double getGlobalTheta(
0075       int endcap, int sector, int subsector, int station, int cscid, bool isME1A, int wiregroup, int halfstrip) const;
0076   double getGlobalThetaFullstrip(
0077       int endcap, int sector, int subsector, int station, int cscid, bool isME1A, int wiregroup, int fullstrip) const;
0078 
0079   // Get sector phi in degrees
0080   double getSectorPhi(int endcap,
0081                       int sector,
0082                       int subsector,
0083                       int station,
0084                       int cscid,
0085                       bool isME1A,
0086                       bool isNeighbor,
0087                       int wiregroup,
0088                       int halfstrip) const;
0089 
0090 private:
0091   const edm::ParameterSet config_;
0092 
0093   int verbose_;
0094 
0095   std::string outdir_;
0096 
0097   bool please_validate_;
0098 
0099   int verbose_sector_;
0100 
0101   bool done_;
0102 
0103   /// Event setup
0104   const CSCGeometry* theCSCGeometry_;
0105 
0106   /// Constants
0107   // [sector_12][station_5][chamber_16]
0108   // NOTE: since Sep 2016, ph_init, ph_cover, ph_disp, th_cover, th_disp are not being used anymore
0109   int ph_init[12][5][16];
0110   int ph_init_full[12][5][16];
0111   int ph_cover[12][5][16];
0112   int ph_disp[12][5][16];
0113   int th_init[12][5][16];
0114   int th_cover[12][5][16];
0115   int th_disp[12][5][16];
0116 
0117   // [station_5][ring_3]
0118   int ph_cover_max[5][3];
0119   int th_cover_max[5][3];
0120 
0121   // [sector_12][station_5][chamber_16][wire_112]
0122   int th_lut[12][5][16][112];
0123   int th_lut_size[12][5][16];
0124 
0125   // [sector_12][station_2][chamber_16][wire_strip_128]  (only ME1/1)
0126   int th_corr_lut[12][2][16][128];
0127   int th_corr_lut_size[12][2][16];
0128 };
0129 
0130 // _____________________________________________________________________________
0131 #define MIN_ENDCAP 1
0132 #define MAX_ENDCAP 2
0133 #define MIN_TRIGSECTOR 1
0134 #define MAX_TRIGSECTOR 6
0135 
0136 #define LOWER_THETA 8.5
0137 #define UPPER_THETA 45.0
0138 
0139 //#define REPRODUCE_OLD_LUTS 1
0140 
0141 MakeCoordLUT::MakeCoordLUT(const edm::ParameterSet& iConfig)
0142     : config_(iConfig),
0143       verbose_(iConfig.getUntrackedParameter<int>("verbosity")),
0144       outdir_(iConfig.getParameter<std::string>("outdir")),
0145       please_validate_(iConfig.getParameter<bool>("please_validate")),
0146       verbose_sector_(2),
0147       done_(false) {
0148   // Zero multi-dimensional arrays
0149   memset(ph_init, 0, sizeof(ph_init));
0150   memset(ph_init_full, 0, sizeof(ph_init_full));
0151   memset(ph_cover, 0, sizeof(ph_cover));
0152   memset(ph_disp, 0, sizeof(ph_disp));
0153   memset(th_init, 0, sizeof(th_init));
0154   memset(th_cover, 0, sizeof(th_cover));
0155   memset(th_disp, 0, sizeof(th_disp));
0156 
0157   memset(ph_cover_max, 0, sizeof(ph_cover_max));
0158   memset(th_cover_max, 0, sizeof(th_cover_max));
0159 
0160   memset(th_lut, 0, sizeof(th_lut));
0161   memset(th_lut_size, 0, sizeof(th_lut_size));
0162 
0163   memset(th_corr_lut, 0, sizeof(th_corr_lut));
0164   memset(th_corr_lut_size, 0, sizeof(th_corr_lut_size));
0165 
0166   assert(CSCConstants::KEY_CLCT_LAYER == CSCConstants::KEY_ALCT_LAYER);
0167 }
0168 
0169 MakeCoordLUT::~MakeCoordLUT() {}
0170 
0171 void MakeCoordLUT::beginRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {
0172   /// Geometry setup
0173   edm::ESHandle<CSCGeometry> cscGeometryHandle;
0174   iSetup.get<MuonGeometryRecord>().get(cscGeometryHandle);
0175   if (!cscGeometryHandle.isValid()) {
0176     std::cout << "ERROR: Unable to get MuonGeometryRecord!" << std::endl;
0177   } else {
0178     theCSCGeometry_ = cscGeometryHandle.product();
0179   }
0180 }
0181 
0182 void MakeCoordLUT::endRun(const edm::Run& iRun, const edm::EventSetup& iSetup) {}
0183 
0184 void MakeCoordLUT::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0185   if (done_)
0186     return;
0187 
0188   generateLUTs();
0189   if (please_validate_)
0190     validateLUTs();
0191   writeFiles();
0192 
0193   done_ = true;
0194   return;
0195 }
0196 
0197 // _____________________________________________________________________________
0198 void MakeCoordLUT::generateLUTs() {
0199   generateLUTs_init();
0200   generateLUTs_run();
0201   generateLUTs_final();
0202 }
0203 
0204 void MakeCoordLUT::generateLUTs_init() {
0205   // Sanity checks
0206   {
0207     // Test ME2/2
0208     CSCDetId id;
0209     id = getCSCDetId(1, 2, 0, 2, 4, false);
0210     assert(id.endcap() == 1 && id.triggerSector() == 2 && id.station() == 2 && id.ring() == 2 &&
0211            CSCTriggerNumbering::triggerCscIdFromLabels(id) == 4);
0212 
0213     // Test ME1/3
0214     id = getCSCDetId(1, 2, 1, 1, 7, false);
0215     assert(id.endcap() == 1 && id.triggerSector() == 2 && id.station() == 1 && id.ring() == 3 &&
0216            CSCTriggerNumbering::triggerCscIdFromLabels(id) == 7);
0217 
0218     // Test ME1/1b
0219     id = getCSCDetId(1, 2, 2, 1, 1, false);
0220     assert(id.endcap() == 1 && id.triggerSector() == 2 && id.station() == 1 && id.ring() == 1 &&
0221            CSCTriggerNumbering::triggerCscIdFromLabels(id) == 1);
0222 
0223     // Test ME1/1a
0224     id = getCSCDetId(1, 2, 2, 1, 1, true);
0225     assert(id.endcap() == 1 && id.triggerSector() == 2 && id.station() == 1 && id.ring() == 4 &&
0226            CSCTriggerNumbering::triggerCscIdFromLabels(id) == 1);
0227   }
0228   return;
0229 }
0230 
0231 // values for ph and th init values hardcoded in verilog zones.v
0232 // these are with offset relative to actual init values to allow for chamber displacement
0233 // [station_5][chamber_16]
0234 // ME1 chambers 13,14,15,16 are neighbor sector chambers 3,6,9,12
0235 // ME2 chambers 10,11 are neighbor sector chambers 3,9
0236 // NOTE: since Sep 2016, th_init_hard and ph_cover_hard are not being used anymore
0237 static const int ph_init_hard[5][16] = {{39, 57, 76, 39, 58, 76, 41, 60, 79, 39, 57, 76, 21, 21, 23, 21},
0238                                         {95, 114, 132, 95, 114, 133, 98, 116, 135, 95, 114, 132, 0, 0, 0, 0},
0239                                         {38, 76, 113, 39, 58, 76, 95, 114, 132, 1, 21, 0, 0, 0, 0, 0},
0240                                         {38, 76, 113, 39, 58, 76, 95, 114, 132, 1, 21, 0, 0, 0, 0, 0},
0241                                         {38, 76, 113, 38, 57, 76, 95, 113, 132, 1, 20, 0, 0, 0, 0, 0}};
0242 
0243 static const int th_init_hard[5][16] = {{1, 1, 1, 42, 42, 42, 94, 94, 94, 1, 1, 1, 1, 42, 94, 1},
0244                                         {1, 1, 1, 42, 42, 42, 94, 94, 94, 1, 1, 1, 0, 0, 0, 0},
0245                                         {1, 1, 1, 48, 48, 48, 48, 48, 48, 1, 48, 0, 0, 0, 0, 0},
0246                                         {1, 1, 1, 40, 40, 40, 40, 40, 40, 1, 40, 0, 0, 0, 0, 0},
0247                                         {2, 2, 2, 34, 34, 34, 34, 34, 34, 2, 34, 0, 0, 0, 0, 0}};
0248 
0249 // hardcoded chamber ph coverage in verilog prim_conv.v
0250 static const int ph_cover_hard[5][16] = {{40, 40, 40, 40, 40, 40, 30, 30, 30, 40, 40, 40, 40, 40, 30, 40},
0251                                          {40, 40, 40, 40, 40, 40, 30, 30, 30, 40, 40, 40, 0, 0, 0, 0},
0252                                          {80, 80, 80, 40, 40, 40, 40, 40, 40, 80, 40, 0, 0, 0, 0, 0},
0253                                          {80, 80, 80, 40, 40, 40, 40, 40, 40, 80, 40, 0, 0, 0, 0, 0},
0254                                          {80, 80, 80, 40, 40, 40, 40, 40, 40, 80, 40, 0, 0, 0, 0, 0}};
0255 
0256 void MakeCoordLUT::generateLUTs_run() {
0257   constexpr double theta_scale = (UPPER_THETA - LOWER_THETA) / 128;  // = 0.28515625 (7 bits encode 128 values)
0258   constexpr double nominal_pitch =
0259       10. / 75.;  // = 0.133333 (ME2/2 strip pitch. 10-degree chamber, 80 strips - 5 overlap strips)
0260 
0261   for (int endcap = MIN_ENDCAP; endcap <= MAX_ENDCAP; ++endcap) {
0262     for (int sector = MIN_TRIGSECTOR; sector <= MAX_TRIGSECTOR; ++sector) {
0263       for (int station = 1; station <= 4; ++station) {
0264         for (int subsector = 0; subsector <= 2; ++subsector) {
0265           for (int chamber = 1; chamber <= 16; ++chamber) {
0266             // ME1 has subsectors 1&2, ME2,3,4 has no subsector (=0)
0267             if ((station == 1 && subsector == 0) || (station != 1 && subsector != 0))
0268               continue;
0269             // Only ME1 subsector 1 has 16 chambers
0270             if (station == 1 && subsector == 2 && chamber > 12)
0271               continue;
0272             // Only ME1 has 12 chambers or more
0273             if (station != 1 && chamber > 11)
0274               continue;
0275 
0276             bool is_me11a = false;
0277             bool is_neighbor = false;
0278 
0279             // Set 'real' CSCID, sector, subsector
0280             int rcscid = chamber;
0281             int rsector = sector;
0282             int rsubsector = subsector;
0283 
0284             if (station == 1) {  // station 1
0285               if (chamber <= 9) {
0286                 rcscid = chamber;
0287               } else if (chamber <= 12) {
0288                 rcscid = (chamber - 9);
0289                 is_me11a = true;
0290               } else if (chamber == 13) {
0291                 rcscid = 3;
0292               } else if (chamber == 14) {
0293                 rcscid = 6;
0294               } else if (chamber == 15) {
0295                 rcscid = 9;
0296               } else if (chamber == 16) {
0297                 rcscid = 3;
0298                 is_me11a = true;
0299               }
0300               if (chamber > 12) {  // is neighbor
0301                 is_neighbor = true;
0302                 rsector = (sector == 1) ? 6 : sector - 1;
0303                 rsubsector = 2;
0304               }
0305 
0306             } else {  // stations 2,3,4
0307               if (chamber <= 9) {
0308                 rcscid = chamber;
0309               } else if (chamber == 10) {
0310                 rcscid = 3;
0311               } else if (chamber == 11) {
0312                 rcscid = 9;
0313               }
0314               if (chamber > 9) {  // is neighbor
0315                 is_neighbor = true;
0316                 rsector = (sector == 1) ? 6 : sector - 1;
0317               }
0318             }
0319 
0320             // Set maxWire, maxStrip
0321             int maxWire = 0;  // can be 32/48/64/96/112
0322             if (station == 1) {
0323               if (rcscid <= 3) {  // ME1/1
0324                 maxWire = 48;
0325               } else if (rcscid <= 6) {  // ME1/2
0326                 maxWire = 64;
0327               } else {  // ME1/3
0328                 maxWire = 32;
0329               }
0330             } else if (station == 2) {
0331               if (rcscid <= 3) {  // ME2/1
0332                 maxWire = 112;
0333               } else {  // ME2/2
0334                 maxWire = 64;
0335               }
0336             } else {              // stations 3,4
0337               if (rcscid <= 3) {  // ME3/1, ME4/1
0338                 maxWire = 96;
0339               } else {  // ME3/2, ME4/2
0340                 maxWire = 64;
0341               }
0342             }
0343 
0344             int maxStrip = 0;  // can be 48/64/80
0345             if (station == 1) {
0346               if (is_me11a) {  // ME1/1a
0347                 maxStrip = 48;
0348               } else if (rcscid <= 3) {  // ME1/1b
0349                 maxStrip = 64;
0350               } else if (6 < rcscid && rcscid <= 9) {  // ME1/3
0351                 maxStrip = 64;
0352               } else {
0353                 maxStrip = 80;
0354               }
0355             } else {
0356               maxStrip = 80;
0357             }
0358 
0359             int topStrip = 0, botStrip = 0, refStrip = 0;
0360             if (station == 1 && rcscid <= 3) {  // ME1/1
0361                                                 // select top and bottom strip according to endcap
0362               // basically, need to hit the corners of the chamber with truncated tilted wires (relevant for ME1/1 only)
0363 #ifdef REPRODUCE_OLD_LUTS
0364               topStrip = (endcap == 2) ? 0 : 47;
0365               botStrip = (endcap == 2) ? 47 : 0;
0366               refStrip = 0;
0367 #else
0368               topStrip = (endcap == 2) ? 0 : maxStrip - 1;
0369               botStrip = (endcap == 2) ? maxStrip - 1 : 0;
0370               refStrip = botStrip;
0371 #endif
0372 
0373             } else {
0374               // take 1/4 of max strip to minimize displacement due to straight wires in polar coordinates (all chambers except ME1/1)
0375               topStrip = maxStrip / 4;
0376               botStrip = maxStrip / 4;
0377               refStrip = botStrip;
0378             }
0379 
0380             const int es = (endcap - 1) * 6 + (sector - 1);
0381             const int st = (station == 1) ? (subsector - 1) : station;
0382             const int ch = (chamber - 1);
0383             assert(es < 12 && st < 5 && ch < 16);
0384             assert(maxWire <= 112);
0385 
0386             // find phi at first and last strips
0387             double fphi_first = getSectorPhi(endcap, rsector, rsubsector, station, rcscid, is_me11a, is_neighbor, 0, 0);
0388             double fphi_last =
0389                 getSectorPhi(endcap, rsector, rsubsector, station, rcscid, is_me11a, is_neighbor, 0, 2 * maxStrip - 1);
0390             double fphi_diff = std::abs(deltaPhiInDegrees(fphi_last, fphi_first)) / 2;  // in double-strip
0391 
0392             // find theta at top and bottom of chamber
0393             double fth_first =
0394                 getGlobalThetaFullstrip(endcap, rsector, rsubsector, station, rcscid, is_me11a, 0, botStrip);
0395             double fth_last =
0396                 getGlobalThetaFullstrip(endcap, rsector, rsubsector, station, rcscid, is_me11a, maxWire - 1, topStrip);
0397 
0398             // find ph_init, ph_init_full, th_init, ph_disp, th_disp constants
0399             int my_ph_init = static_cast<int>(std::round(fphi_first / nominal_pitch));
0400             int my_ph_init_full = static_cast<int>(std::round(fphi_first / (nominal_pitch / 8.)));  // 1/8-strip pitch
0401             int my_ph_cover = static_cast<int>(std::round(fphi_diff / nominal_pitch));
0402             int my_th_init = static_cast<int>(std::round((fth_first - LOWER_THETA) / theta_scale));
0403             int my_th_cover = static_cast<int>(std::round((fth_last - fth_first) / theta_scale));
0404 
0405             // calculate displacements from hardcoded init values
0406             int my_ph_disp = (my_ph_init / 2 - 2 * ph_init_hard[st][ch]);  // in double-strip
0407             if (deltaPhiInDegrees(fphi_first, fphi_last) > 0.)
0408               my_ph_disp -= ph_cover_hard[st][ch];
0409             int my_th_disp = (my_th_init - th_init_hard[st][ch]);
0410 
0411 #ifdef REPRODUCE_OLD_LUTS
0412             // widen ME1/1 coverage slightly, because of odd geometry of truncated wiregroups
0413             if (station == 1 && rcscid <= 3) {  // ME1/1
0414               my_th_cover += 2;
0415             }
0416 #endif
0417 
0418             ph_init[es][st][ch] = my_ph_init;
0419             ph_init_full[es][st][ch] = my_ph_init_full;
0420             ph_cover[es][st][ch] = my_ph_cover;
0421             ph_disp[es][st][ch] = my_ph_disp;
0422             th_init[es][st][ch] = my_th_init;
0423             th_cover[es][st][ch] = my_th_cover;
0424             th_disp[es][st][ch] = my_th_disp;
0425 
0426             if (verbose_ > 0 && sector == verbose_sector_) {
0427               double fphi_first_global = getGlobalPhi(endcap, rsector, rsubsector, station, rcscid, is_me11a, 0, 0);
0428               std::cout << "::generateLUTs_run()"
0429                         << " -- endcap " << endcap << " sec " << sector << " st " << st << " ch " << ch + 1
0430                         << " maxWire " << maxWire << " maxStrip " << maxStrip
0431                         << " -- fphi_first_global: " << fphi_first_global << " fphi_first: " << fphi_first
0432                         << " fphi_last: " << fphi_last << " fth_first: " << fth_first << " fth_last: " << fth_last
0433                         << " ph_init: " << my_ph_init << " ph_init_full: " << my_ph_init_full
0434                         << " ph_cover: " << my_ph_cover << " ph_disp: " << my_ph_disp << " th_init: " << my_th_init
0435                         << " th_cover: " << my_th_cover << " th_disp: " << my_th_disp << std::endl;
0436             }
0437 
0438             // make LUT for wire -> theta
0439             for (int wire = 0; wire < maxWire; ++wire) {
0440               double fth_wire =
0441                   getGlobalThetaFullstrip(endcap, rsector, rsubsector, station, rcscid, is_me11a, wire, refStrip);
0442               double fth_diff = fth_wire - fth_first;
0443               int th_diff = static_cast<int>(std::round(fth_diff / theta_scale));
0444               assert(th_diff >= 0);
0445 
0446               if (wire == 0)
0447                 th_lut_size[es][st][ch] = maxWire;
0448               th_lut[es][st][ch][wire] = th_diff;
0449 
0450               if (verbose_ > 0 && sector == verbose_sector_) {
0451                 std::cout << "::generateLUTs_run()"
0452                           << " -- endcap " << endcap << " sec " << sector << " st " << st << " ch " << ch + 1
0453                           << " wire " << wire << " strip " << refStrip << " -- fth_first: " << fth_first
0454                           << " fth_wire: " << fth_wire << " fth_diff: " << fth_diff << " th_diff: " << th_diff
0455                           << std::endl;
0456               }
0457             }  // end loop over wire
0458 
0459             // make LUT for (wire,strip) index -> theta correction for ME1/1 where the wires are tilted
0460             if (station == 1 && rcscid <= 3 && !is_me11a) {  // ME1/1b
0461               assert(maxWire == 48 && maxStrip == 64);       // ME1/1b
0462 
0463               // select correction points at 1/6, 3/6 and 5/6 of chamber wg range
0464               // this makes construction of LUT address in firmware much easier
0465               int index = 0;
0466 
0467               for (int wire = maxWire / 6; wire < maxWire; wire += maxWire / 3) {
0468                 double fth0 =
0469                     getGlobalThetaFullstrip(endcap, rsector, rsubsector, station, rcscid, is_me11a, wire, refStrip);
0470 
0471                 // pattern search works in double-strip, so take every other strip
0472                 for (int strip = 0; strip < maxStrip; strip += 2) {
0473                   double fth1 =
0474                       getGlobalThetaFullstrip(endcap, rsector, rsubsector, station, rcscid, is_me11a, wire, strip);
0475                   double fth_diff = fth1 - fth0;
0476 
0477 #ifdef REPRODUCE_OLD_LUTS
0478                   // for chambers in negative endcap, the wire tilt is the opposite way
0479                   fth_diff = (endcap == 2) ? -fth_diff : fth_diff;
0480 #endif
0481 
0482                   int th_diff = static_cast<int>(std::round(fth_diff / theta_scale));
0483                   assert(th_diff >= 0);
0484                   assert(index <= 96);  // (3) [wire] x (64/2) [strip]
0485 
0486                   if (index == 0)
0487                     th_corr_lut_size[es][st][ch] = 96;  // (3) [wire] x (64/2) [strip]
0488                   th_corr_lut[es][st][ch][index] = th_diff;
0489 
0490                   if (verbose_ > 0 && sector == verbose_sector_) {
0491                     std::cout << "::generateLUTs_run()"
0492                               << " -- endcap " << endcap << " sec " << sector << " st " << st << " ch " << ch + 1
0493                               << " wire " << wire << " strip " << strip << " -- fth0: " << fth0 << " fth1: " << fth1
0494                               << " fth_diff: " << fth_diff << " th_diff: " << th_diff << std::endl;
0495                   }
0496 
0497                   ++index;
0498                 }  // end loop over strip
0499               }    // end loop over wire
0500             }      // end if ME1/1b
0501           }        // end loop over chamber
0502         }          // end loop over subsector
0503       }            // end loop over station
0504     }              // end loop over sector
0505   }                // end loop over endcap
0506   return;
0507 }
0508 
0509 void MakeCoordLUT::generateLUTs_final() {
0510   // update max coverages
0511   for (int es = 0; es < 12; ++es) {
0512     for (int st = 0; st < 5; ++st) {
0513       for (int ch = 0; ch < 16; ++ch) {
0514         if (ch > 9 - 1)
0515           continue;  // exclude neighbors, exclude ME1/1a
0516 
0517         int ch_type = ch / 3;
0518         if (st > 1 && ch_type > 1)
0519           ch_type = 1;  // stations 2,3,4 have only 2 chamber types (a.k.a rings)
0520 
0521         if (ph_cover_max[st][ch_type] < ph_cover[es][st][ch])
0522           ph_cover_max[st][ch_type] = ph_cover[es][st][ch];
0523         if (th_cover_max[st][ch_type] < th_cover[es][st][ch])
0524           th_cover_max[st][ch_type] = th_cover[es][st][ch];
0525       }  // end loop over ch
0526     }    // end loop over st
0527   }      // end loop over es
0528 
0529   for (int st = 0; st < 5; ++st) {
0530     for (int ch_type = 0; ch_type < 3; ++ch_type) {
0531       if (verbose_ > 0) {
0532         std::cout << "::generateLUTs_final()"
0533                   << " -- st " << st << " ch_type " << ch_type + 1 << " -- ph_cover_max: " << ph_cover_max[st][ch_type]
0534                   << " th_cover_max: " << th_cover_max[st][ch_type] << std::endl;
0535       }
0536     }  // end loop over ch_type
0537   }    // end loop over st
0538   return;
0539 }
0540 
0541 // Compare simulated (with floating-point) vs emulated (fixed-point) phi and theta coordinates
0542 void MakeCoordLUT::validateLUTs() {
0543   std::stringstream filename;
0544 
0545   filename << outdir_ << "/"
0546            << "validate.root";
0547   TFile* tfile = TFile::Open(filename.str().c_str(), "RECREATE");
0548   filename.str("");
0549   filename.clear();
0550 
0551   // Create TTree
0552   int lut_id = 0;
0553   int es = 0;
0554   int st = 0;
0555   int ch = 0;
0556   //
0557   int endcap = 0;
0558   int station = 0;
0559   int sector = 0;
0560   int subsector = 0;
0561   int ring = 0;
0562   int chamber = 0;
0563   int CSC_ID = 0;
0564   //
0565   int strip = 0;  // it is half-strip, despite the name
0566   int wire = 0;   // it is wiregroup, despite the name
0567   int fph_int = 0;
0568   int fth_int = 0;
0569   double fph_emu = 0.;  // in degrees
0570   double fth_emu = 0.;  // in degrees
0571   double fph_sim = 0.;  // in degrees
0572   double fth_sim = 0.;  // in degrees
0573 
0574   TTree* ttree = new TTree("tree", "tree");
0575   ttree->Branch("lut_id", &lut_id);
0576   ttree->Branch("es", &es);
0577   ttree->Branch("st", &st);
0578   ttree->Branch("ch", &ch);
0579   //
0580   ttree->Branch("endcap", &endcap);
0581   ttree->Branch("station", &station);
0582   ttree->Branch("sector", &sector);
0583   ttree->Branch("subsector", &subsector);
0584   ttree->Branch("ring", &ring);
0585   ttree->Branch("chamber", &chamber);
0586   ttree->Branch("CSC_ID", &CSC_ID);
0587   //
0588   ttree->Branch("strip", &strip);
0589   ttree->Branch("wire", &wire);
0590   ttree->Branch("fph_int", &fph_int);
0591   ttree->Branch("fth_int", &fth_int);
0592   ttree->Branch("fph_emu", &fph_emu);
0593   ttree->Branch("fth_emu", &fth_emu);
0594   ttree->Branch("fph_sim", &fph_sim);
0595   ttree->Branch("fth_sim", &fth_sim);
0596 
0597   for (es = 0; es < 12; ++es) {
0598     for (lut_id = 0; lut_id < 61; ++lut_id) {  // every sector has 61 LUT id
0599       // Retrieve st, ch from lut_id
0600       if (lut_id < 16) {
0601         st = 0;
0602         ch = lut_id - 0;
0603       } else if (lut_id < 28) {
0604         st = 1;
0605         ch = lut_id - 16;
0606       } else if (lut_id < 39) {
0607         st = 2;
0608         ch = lut_id - 28;
0609       } else if (lut_id < 50) {
0610         st = 3;
0611         ch = lut_id - 39;
0612       } else {
0613         st = 4;
0614         ch = lut_id - 50;
0615       }
0616       assert(es >= 0 && st >= 0 && ch >= 0);
0617       assert(es < 12 && st < 5 && ch < 16);
0618 
0619       // Retrieve endcap, sector, subsector, station, chamber
0620       endcap = (es / 6) + 1;
0621       sector = (es % 6) + 1;
0622       subsector = (st <= 1) ? st + 1 : 0;
0623       station = (st <= 1) ? 1 : st;
0624       chamber = ch + 1;
0625 
0626       bool is_me11a = false;
0627       bool is_neighbor = false;
0628 
0629       // Set 'real' CSCID, sector, subsector
0630       int rcscid = chamber;
0631       int rsector = sector;
0632       int rsubsector = subsector;
0633 
0634       if (station == 1) {  // station 1
0635         if (chamber <= 9) {
0636           rcscid = chamber;
0637         } else if (chamber <= 12) {
0638           rcscid = (chamber - 9);
0639           is_me11a = true;
0640         } else if (chamber == 13) {
0641           rcscid = 3;
0642         } else if (chamber == 14) {
0643           rcscid = 6;
0644         } else if (chamber == 15) {
0645           rcscid = 9;
0646         } else if (chamber == 16) {
0647           rcscid = 3;
0648           is_me11a = true;
0649         }
0650         if (chamber > 12) {  // is neighbor
0651           is_neighbor = true;
0652           rsector = (sector == 1) ? 6 : sector - 1;
0653           rsubsector = 2;
0654         }
0655 
0656       } else {  // stations 2,3,4
0657         if (chamber <= 9) {
0658           rcscid = chamber;
0659         } else if (chamber == 10) {
0660           rcscid = 3;
0661         } else if (chamber == 11) {
0662           rcscid = 9;
0663         }
0664         if (chamber > 9) {  // is neighbor
0665           is_neighbor = true;
0666           rsector = (sector == 1) ? 6 : sector - 1;
0667         }
0668       }
0669 
0670       CSC_ID = rcscid;
0671 
0672       // Set maxWire, maxStrip
0673       const CSCDetId cscDetId = getCSCDetId(endcap, rsector, rsubsector, station, rcscid, is_me11a);
0674       const CSCChamber* chamb = theCSCGeometry_->chamber(cscDetId);
0675       const CSCLayerGeometry* layerGeom = chamb->layer(CSCConstants::KEY_CLCT_LAYER)->geometry();
0676 
0677       ring = cscDetId.ring();
0678       const int maxWire = layerGeom->numberOfWireGroups();
0679       const int maxStrip = layerGeom->numberOfStrips();
0680 
0681       // _______________________________________________________________________
0682       // Copied from PrimitiveConversion
0683 
0684       const int fw_endcap = (es / 6);
0685       const int fw_sector = (es % 6);
0686       const int fw_station = st;
0687       const int fw_cscid = is_me11a ? (is_neighbor ? ch - 3 : ch - 9) : ch;
0688 
0689       // Is this chamber mounted in reverse direction?
0690       bool ph_reverse = false;
0691       if ((fw_endcap == 0 && fw_station >= 3) || (fw_endcap == 1 && fw_station < 3))
0692         ph_reverse = true;
0693 
0694       // Is this 10-deg or 20-deg chamber?
0695       bool is_10degree = false;
0696       if ((fw_station <= 1) ||                                                       // ME1
0697           (fw_station >= 2 && ((fw_cscid >= 3 && fw_cscid <= 8) || fw_cscid == 10))  // ME2,3,4/2
0698       ) {
0699         is_10degree = true;
0700       }
0701 
0702       assert(ph_reverse ==
0703              isStripPhiCounterClockwise(getCSCDetId(endcap, rsector, rsubsector, station, rcscid, is_me11a)));
0704 
0705       for (wire = 0; wire < maxWire; ++wire) {
0706         for (strip = 0; strip < 2 * maxStrip; ++strip) {
0707           const int fw_strip = strip;  // it is half-strip, despite the name
0708           const int fw_wire = wire;    // it is wiregroup, despite the name
0709 
0710           // ___________________________________________________________________
0711           // phi conversion
0712 
0713           // Convert half-strip into 1/8-strip
0714           int eighth_strip = 0;
0715 
0716           // Apply phi correction from CLCT pattern number
0717           int clct_pat_corr = 0;
0718           int clct_pat_corr_sign = 1;
0719 
0720           if (is_10degree) {
0721             eighth_strip = fw_strip << 2;  // full precision, uses only 2 bits of pattern correction
0722             eighth_strip += clct_pat_corr_sign * (clct_pat_corr >> 1);
0723           } else {
0724             eighth_strip = fw_strip << 3;  // multiply by 2, uses all 3 bits of pattern correction
0725             eighth_strip += clct_pat_corr_sign * (clct_pat_corr >> 0);
0726           }
0727 
0728           // Multiplicative factor for eighth_strip
0729           int factor = 1024;
0730           if (station == 1 && ring == 4)
0731             factor = 1707;  // ME1/1a
0732           else if (station == 1 && ring == 1)
0733             factor = 1301;  // ME1/1b
0734           else if (station == 1 && ring == 3)
0735             factor = 947;  // ME1/3
0736 
0737           // ph_tmp is full-precision phi, but local to chamber (counted from strip 0)
0738           // full phi precision: 0.016666 deg (1/8-strip)
0739           // zone phi precision: 0.533333 deg (4-strip, 32 times coarser than full phi precision)
0740           int ph_tmp = (eighth_strip * factor) >> 10;
0741           int ph_tmp_sign = (ph_reverse == 0) ? 1 : -1;
0742 
0743           int fph = ph_init_full[es][st][ch];
0744           fph = fph + ph_tmp_sign * ph_tmp;
0745 
0746           // ph_init_hard is used to calculate zone_hit in the firmware
0747           assert(((fph + (1 << 4)) >> 5) >= ph_init_hard[st][ch]);
0748 
0749           // ___________________________________________________________________
0750           // theta conversion
0751 
0752           // Make ME1/1a the same as ME1/1b when using th_lut and th_corr_lut
0753           int ch2 = is_me11a ? (is_neighbor ? ch - 3 : ch - 9) : ch;
0754 
0755           // th_tmp is theta local to chamber
0756           int pc_wire_id = (fw_wire & 0x7f);  // 7-bit
0757           assert(pc_wire_id < th_lut_size[es][st][ch2]);
0758           int th_tmp = th_lut[es][st][ch2][pc_wire_id];
0759 
0760           // For ME1/1 with tilted wires, add theta correction as a function of (wire,strip) index
0761           if (station == 1 && (ring == 1 || ring == 4)) {
0762 #ifdef REPRODUCE_OLD_LUTS
0763             int pc_wire_strip_id =
0764                 (((fw_wire >> 4) & 0x3) << 5) | ((eighth_strip >> 4) & 0x1f);  // 2-bit from wire, 5-bit from 2-strip
0765             assert(pc_wire_strip_id < th_corr_lut_size[es][st][ch2]);
0766             int th_corr = th_corr_lut[es][st][ch2][pc_wire_strip_id];
0767             int th_corr_sign = (ph_reverse == 0) ? 1 : -1;
0768 
0769             th_tmp = th_tmp + th_corr_sign * th_corr;
0770 
0771             // Check that correction did not make invalid value outside chamber coverage
0772             const int th_negative = 50;
0773             const int th_coverage = 45;
0774 
0775             if (th_tmp > th_negative || th_tmp < 0 || fw_wire == 0)
0776               th_tmp = 0;  // limit at the bottom
0777             if (th_tmp > th_coverage)
0778               th_tmp = th_coverage;  // limit at the top
0779 #else
0780             int pc_wire_strip_id =
0781                 (((fw_wire >> 4) & 0x3) << 5) | ((eighth_strip >> 4) & 0x1f);  // 2-bit from wire, 5-bit from 2-strip
0782             if (is_me11a)
0783               pc_wire_strip_id =
0784                   (((fw_wire >> 4) & 0x3) << 5) |
0785                   ((((eighth_strip * 341) >> 8) >> 4) & 0x1f);  // correct for ME1/1a strip number (341/256 =~ 1.333)
0786             assert(pc_wire_strip_id < th_corr_lut_size[es][st][ch2]);
0787             int th_corr = th_corr_lut[es][st][ch2][pc_wire_strip_id];
0788 
0789             th_tmp = th_tmp + th_corr;
0790             assert(th_tmp >= 0);
0791 
0792             // Check that correction did not make invalid value outside chamber coverage
0793             const int th_coverage = 46;  // max coverage for front chamber is 47, max coverage for rear chamber is 45
0794 
0795             if (fw_wire == 0)
0796               th_tmp = 0;  // limit at the bottom
0797             if (th_tmp > th_coverage)
0798               th_tmp = th_coverage;  // limit at the top
0799 #endif
0800           }
0801 
0802           // theta precision: 0.28515625 deg
0803           int th = th_init[es][st][ch];
0804           th = th + th_tmp;
0805 
0806           // Protect against invalid value
0807           th = (th == 0) ? 1 : th;
0808 
0809           // ___________________________________________________________________
0810           // Finally
0811 
0812           // emulated phi and theta coordinates from fixed-point operations
0813           fph_int = fph;
0814           fph_emu = static_cast<double>(fph_int);
0815           fph_emu = fph_emu / 60.;
0816           fph_emu = fph_emu - 22. + 15. + (60. * fw_sector);
0817           fph_emu = deltaPhiInDegrees(fph_emu, 0.);  // reduce to [-180,180]
0818 
0819           fth_int = th;
0820           fth_emu = static_cast<double>(fth_int);
0821           fth_emu = (fth_emu * (45.0 - 8.5) / 128. + 8.5);
0822 
0823           // simulated phi and theta coordinates from floating-point operations
0824           fph_sim = getGlobalPhi(endcap, rsector, rsubsector, station, rcscid, is_me11a, wire, strip);
0825           fth_sim = getGlobalTheta(endcap, rsector, rsubsector, station, rcscid, is_me11a, wire, strip);
0826 
0827           ttree->Fill();
0828 
0829           if (verbose_ > 1 && sector == verbose_sector_) {
0830             std::cout << "::validateLUTs()"
0831                       << " -- endcap " << endcap << " sec " << sector << " st " << st << " ch " << ch + 1 << " wire "
0832                       << wire << " strip " << strip << " -- fph_int: " << fph_int << " fph_emu: " << fph_emu
0833                       << " fph_sim: " << fph_sim << " -- fth_int: " << fth_int << " fth_emu: " << fth_emu
0834                       << " fth_sim: " << fth_sim << std::endl;
0835           }
0836         }  // end loop over strip
0837       }    // end loop over wire
0838     }      // end loop over lut_id
0839   }        // end loop over es
0840 
0841   ttree->Write();
0842   tfile->Close();
0843   return;
0844 }
0845 
0846 // produce the LUT text files
0847 void MakeCoordLUT::writeFiles() {
0848   int num_of_files = 0;
0849 
0850   std::stringstream filename;
0851 
0852   for (int es = 0; es < 12; ++es) {
0853     int endcap = (es / 6) + 1;
0854     int sector = (es % 6) + 1;
0855 
0856     // write files: ph_init, ph_init_full, th_init, ph_disp, th_disp
0857     std::ofstream ph_init_fs;
0858     filename << outdir_ << "/"
0859              << "ph_init_endcap_" << endcap << "_sect_" << sector << ".lut";
0860     ph_init_fs.open(filename.str().c_str());
0861     filename.str("");
0862     filename.clear();
0863 
0864     std::ofstream th_init_fs;
0865     filename << outdir_ << "/"
0866              << "th_init_endcap_" << endcap << "_sect_" << sector << ".lut";
0867     th_init_fs.open(filename.str().c_str());
0868     filename.str("");
0869     filename.clear();
0870 
0871     std::ofstream ph_disp_fs;
0872     filename << outdir_ << "/"
0873              << "ph_disp_endcap_" << endcap << "_sect_" << sector << ".lut";
0874     ph_disp_fs.open(filename.str().c_str());
0875     filename.str("");
0876     filename.clear();
0877 
0878     std::ofstream th_disp_fs;
0879     filename << outdir_ << "/"
0880              << "th_disp_endcap_" << endcap << "_sect_" << sector << ".lut";
0881     th_disp_fs.open(filename.str().c_str());
0882     filename.str("");
0883     filename.clear();
0884 
0885     for (int st = 0; st < 5; ++st) {
0886       const int max_ch = (st == 0) ? 16 : (st == 1) ? 12 : 11;
0887 
0888       std::ofstream ph_init_full_fs;
0889       filename << outdir_ << "/"
0890                << "ph_init_full_endcap_" << endcap << "_sect_" << sector << "_st_" << st << ".lut";
0891       ph_init_full_fs.open(filename.str().c_str());
0892       filename.str("");
0893       filename.clear();
0894 
0895       for (int ch = 0; ch < max_ch; ++ch) {
0896         assert(es < 12 && st < 5 && ch < 16);
0897 
0898         ph_init_fs << std::hex << ph_init[es][st][ch] << std::endl;
0899         ph_init_full_fs << std::hex << ph_init_full[es][st][ch] << std::endl;
0900         th_init_fs << std::hex << th_init[es][st][ch] << std::endl;
0901         ph_disp_fs << std::hex << ph_disp[es][st][ch] << std::endl;
0902         th_disp_fs << std::hex << th_disp[es][st][ch] << std::endl;
0903       }  // end loop over ch
0904 
0905       ph_init_full_fs.close();
0906       ++num_of_files;
0907     }  // end loop over st
0908 
0909     ph_init_fs.close();
0910     ++num_of_files;
0911     th_init_fs.close();
0912     ++num_of_files;
0913     ph_disp_fs.close();
0914     ++num_of_files;
0915     th_disp_fs.close();
0916     ++num_of_files;
0917 
0918     // write files: th_lut, th_corr_lut
0919     for (int st = 0; st < 5; ++st) {
0920       const int max_ch = (st == 0) ? 16 : (st == 1) ? 12 : 11;
0921 
0922       for (int ch = 0; ch < max_ch; ++ch) {
0923         assert(es < 12 && st < 5 && ch < 16);
0924 
0925         int subsector = (st <= 1) ? st + 1 : 0;
0926         int station = (st <= 1) ? 1 : st;
0927         int chamber = ch + 1;
0928 
0929         std::ofstream th_lut_fs;
0930         if (station == 1) {
0931           filename << outdir_ << "/"
0932                    << "vl_th_lut_endcap_" << endcap << "_sec_" << sector << "_sub_" << subsector << "_st_" << station
0933                    << "_ch_" << chamber << ".lut";
0934         } else {
0935           filename << outdir_ << "/"
0936                    << "vl_th_lut_endcap_" << endcap << "_sec_" << sector << "_st_" << station << "_ch_" << chamber
0937                    << ".lut";
0938         }
0939         th_lut_fs.open(filename.str().c_str());
0940         filename.str("");
0941         filename.clear();
0942 
0943         const int maxWire = th_lut_size[es][st][ch];
0944         for (int wire = 0; wire < maxWire; ++wire) {
0945           th_lut_fs << std::hex << th_lut[es][st][ch][wire] << std::endl;
0946         }
0947         th_lut_fs.close();
0948         ++num_of_files;
0949 
0950         if (station == 1 && (ch == 0 || ch == 1 || ch == 2 || ch == 12)) {  // ME1/1 chambers
0951           std::ofstream th_corr_lut_fs;
0952           filename << outdir_ << "/"
0953                    << "vl_th_corr_lut_endcap_" << endcap << "_sec_" << sector << "_sub_" << subsector << "_st_"
0954                    << station << "_ch_" << ch + 1 << ".lut";
0955           th_corr_lut_fs.open(filename.str().c_str());
0956           filename.str("");
0957           filename.clear();
0958 
0959           const int n = th_corr_lut_size[es][st][ch];
0960           for (int index = 0; index < n; ++index) {
0961             th_corr_lut_fs << std::hex << th_corr_lut[es][st][ch][index] << std::endl;
0962           }
0963           th_corr_lut_fs.close();
0964           ++num_of_files;
0965         }
0966       }  // end loop over ch
0967     }    // end loop over st
0968 
0969   }  // end loop over es
0970 
0971   std::cout << "[INFO] Generated " << num_of_files << " LUT files." << std::endl;
0972 
0973   // Expect 12 sectors x (7 th_corr_lut + 61 th_lut + 4 ph_init/th_init/ph_disp/th_disp + 5 ph_init_full)
0974   assert(num_of_files == 12 * (7 + 61 + 4 + 5));
0975   return;
0976 }
0977 
0978 // _____________________________________________________________________________
0979 CSCDetId MakeCoordLUT::getCSCDetId(int endcap, int sector, int subsector, int station, int cscid, bool isME1A) const {
0980   int ring = isME1A ? 4 : CSCTriggerNumbering::ringFromTriggerLabels(station, cscid);
0981   int chamber = CSCTriggerNumbering::chamberFromTriggerLabels(sector, subsector, station, cscid);
0982   const CSCDetId cscDetId = CSCDetId(endcap, station, ring, chamber, CSCConstants::KEY_CLCT_LAYER);
0983   return cscDetId;
0984 }
0985 
0986 bool MakeCoordLUT::isStripPhiCounterClockwise(const CSCDetId& cscDetId) const {
0987   const CSCChamber* chamb = theCSCGeometry_->chamber(cscDetId);
0988   const CSCLayer* layer = chamb->layer(CSCConstants::KEY_CLCT_LAYER);
0989 
0990   const double phi1 = layer->centerOfStrip(1).phi();
0991   const double phi2 = layer->centerOfStrip(2).phi();
0992   bool ccw = (deltaPhiInRadians(phi1, phi2) > 0.);
0993   return ccw;
0994 }
0995 
0996 double MakeCoordLUT::getStripPitch(const CSCDetId& cscDetId) const {
0997   const CSCChamber* chamb = theCSCGeometry_->chamber(cscDetId);
0998   const CSCLayerGeometry* layerGeom = chamb->layer(CSCConstants::KEY_CLCT_LAYER)->geometry();
0999 
1000   constexpr double _rad_to_deg = 180. / M_PI;
1001   double pitch = layerGeom->stripPhiPitch() * _rad_to_deg;
1002   return pitch;
1003 }
1004 
1005 double MakeCoordLUT::getGlobalPhi(
1006     int endcap, int sector, int subsector, int station, int cscid, bool isME1A, int wiregroup, int halfstrip) const {
1007   int fullstrip = (halfstrip / 2);
1008   int oddhs = (halfstrip % 2);
1009   double phi = getGlobalPhiFullstrip(endcap, sector, subsector, station, cscid, isME1A, wiregroup, fullstrip);
1010 
1011   // Add half-strip offset
1012   // strip width/4 gives the offset of the half-strip center w.r.t the strip center
1013   const CSCDetId cscDetId = getCSCDetId(endcap, sector, subsector, station, cscid, isME1A);
1014   double pitch = getStripPitch(cscDetId) / 4.0;
1015   bool ph_reverse = isStripPhiCounterClockwise(cscDetId);
1016 
1017   pitch = (ph_reverse == 1) ? -pitch : pitch;  // subtract half-strip if phi decreases as strip number increases
1018   pitch = (oddhs == 0) ? -pitch : pitch;       // subtract even half-strip or add odd half-strip
1019   phi += pitch;
1020   return phi;
1021 }
1022 
1023 double MakeCoordLUT::getGlobalPhiFullstrip(
1024     int endcap, int sector, int subsector, int station, int cscid, bool isME1A, int wiregroup, int fullstrip) const {
1025   const CSCDetId cscDetId = getCSCDetId(endcap, sector, subsector, station, cscid, isME1A);
1026   const CSCChamber* chamb = theCSCGeometry_->chamber(cscDetId);
1027   const CSCLayerGeometry* layerGeom = chamb->layer(CSCConstants::KEY_CLCT_LAYER)->geometry();
1028 
1029   const LocalPoint& lp = layerGeom->stripWireGroupIntersection(
1030       fullstrip + 1, wiregroup + 1);  // strip and wg in geometry routines start from 1
1031   const GlobalPoint& gp = chamb->layer(CSCConstants::KEY_ALCT_LAYER)->surface().toGlobal(lp);
1032 
1033   constexpr double _rad_to_deg = 180. / M_PI;
1034   double phi = gp.barePhi() * _rad_to_deg;
1035   return phi;
1036 }
1037 
1038 double MakeCoordLUT::getGlobalTheta(
1039     int endcap, int sector, int subsector, int station, int cscid, bool isME1A, int wiregroup, int halfstrip) const {
1040   int fullstrip = (halfstrip / 2);
1041   return getGlobalThetaFullstrip(endcap, sector, subsector, station, cscid, isME1A, wiregroup, fullstrip);
1042 }
1043 
1044 double MakeCoordLUT::getGlobalThetaFullstrip(
1045     int endcap, int sector, int subsector, int station, int cscid, bool isME1A, int wiregroup, int fullstrip) const {
1046   const CSCDetId cscDetId = getCSCDetId(endcap, sector, subsector, station, cscid, isME1A);
1047   const CSCChamber* chamb = theCSCGeometry_->chamber(cscDetId);
1048   const CSCLayerGeometry* layerGeom = chamb->layer(CSCConstants::KEY_CLCT_LAYER)->geometry();
1049 
1050   //const LocalPoint& lp  = layerGeom->stripWireGroupIntersection(fullstrip+1, wiregroup+1); // strip and wg in geometry routines start from 1
1051   const LocalPoint& lp =
1052       layerGeom->intersectionOfStripAndWire(fullstrip + 1, layerGeom->middleWireOfGroup(wiregroup + 1));
1053   const GlobalPoint& gp = chamb->layer(CSCConstants::KEY_ALCT_LAYER)->surface().toGlobal(lp);
1054 
1055   constexpr double _rad_to_deg = 180. / M_PI;
1056   double theta = gp.theta() * _rad_to_deg;
1057   theta = (endcap == 2) ? (180. - theta) : theta;  // put theta in the range of 0 - 180 degrees for negative endcap
1058   return theta;
1059 }
1060 
1061 double MakeCoordLUT::getSectorPhi(int endcap,
1062                                   int sector,
1063                                   int subsector,
1064                                   int station,
1065                                   int cscid,
1066                                   bool isME1A,
1067                                   bool isNeighbor,
1068                                   int wiregroup,
1069                                   int halfstrip) const {
1070   double globalPhi = getGlobalPhi(endcap, sector, subsector, station, cscid, isME1A, wiregroup, halfstrip);
1071 
1072   // sector boundary should not depend on station, cscid, etc. For now, take station 2 csc 1 strip 0 as boundary, -2 deg (Darin, 2009-09-18)
1073   // correction for sector overlap: take sector boundary at neighbor sector station 2 csc 3 strip 0, - 2 deg (Matt, 2016-03-07)
1074   int sector_n = (sector == 1) ? 6 : sector - 1;  // neighbor sector
1075   if (isNeighbor) {
1076     sector_n = sector;  // same sector
1077   }
1078 
1079   const int maxStrip = 80;  // ME2
1080   const int firstWire = 0;
1081   const int firstStrip = (endcap == 1) ? 0 : maxStrip - 1;
1082   double sectorStartPhi = getGlobalPhiFullstrip(endcap, sector_n, 0, 2, 3, false, firstWire, firstStrip) - 2.;
1083 
1084 #ifndef REPRODUCE_OLD_LUTS
1085   // but sector boundary does depend on endcap. apply additional correction to make integer phi 0
1086   // lines up at -22 deg (Jia Fu, 2016-11-12)
1087   //sectorStartPhi = (endcap == 2) ? sectorStartPhi + 36./60 : sectorStartPhi + 28./60;
1088 
1089   // Manually lines up at -22 deg (Jia Fu, 2018-09-19)
1090   double oldSectorStartPhi = sectorStartPhi;
1091   sectorStartPhi = -22. + 15. + (60. * (sector - 1));
1092   if (isNeighbor) {
1093     // This chamber comes from the neighbor sector into the native sector
1094     // Use the native sector sectorStartPhi (+60 deg)
1095     sectorStartPhi += 60.;
1096   }
1097   if (sectorStartPhi > 180.)
1098     sectorStartPhi -= 360.;
1099   assert(std::abs(oldSectorStartPhi - sectorStartPhi) < 2.);  // sanity check
1100 #endif
1101 
1102   double res = deltaPhiInDegrees(globalPhi, sectorStartPhi);
1103   assert(res >= 0.);
1104   return res;
1105 }
1106 
1107 // DEFINE THIS AS A PLUG-IN
1108 #include "FWCore/Framework/interface/MakerMacros.h"
1109 DEFINE_FWK_MODULE(MakeCoordLUT);