Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-09-07 04:36:59 

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

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!