Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:19:48

0001 #include "L1Trigger/DTTriggerPhase2/interface/MuonPathConfirmator.h"
0002 #include <cmath>
0003 #include <memory>
0004 
0005 using namespace edm;
0006 using namespace std;
0007 using namespace cmsdt;
0008 
0009 // ============================================================================
0010 // Constructors and destructor
0011 // ============================================================================
0012 MuonPathConfirmator::MuonPathConfirmator(const ParameterSet &pset, edm::ConsumesCollector &iC)
0013     : debug_(pset.getUntrackedParameter<bool>("debug")),
0014       minx_match_2digis_(pset.getParameter<double>("minx_match_2digis")) {
0015   if (debug_)
0016     LogDebug("MuonPathConfirmator") << "MuonPathConfirmator: constructor";
0017 
0018   //shift phi
0019   int rawId;
0020   shift_filename_ = pset.getParameter<edm::FileInPath>("shift_filename");
0021   std::ifstream ifin3(shift_filename_.fullPath());
0022   double shift;
0023   if (ifin3.fail()) {
0024     throw cms::Exception("Missing Input File")
0025         << "MuonPathConfirmator::MuonPathConfirmator() -  Cannot find " << shift_filename_.fullPath();
0026   }
0027   while (ifin3.good()) {
0028     ifin3 >> rawId >> shift;
0029     shiftinfo_[rawId] = shift;
0030   }
0031 
0032   int wh, st, se, maxdrift;
0033   maxdrift_filename_ = pset.getParameter<edm::FileInPath>("maxdrift_filename");
0034   std::ifstream ifind(maxdrift_filename_.fullPath());
0035   if (ifind.fail()) {
0036     throw cms::Exception("Missing Input File")
0037         << "MPSLFilter::MPSLFilter() -  Cannot find " << maxdrift_filename_.fullPath();
0038   }
0039   while (ifind.good()) {
0040     ifind >> wh >> st >> se >> maxdrift;
0041     maxdriftinfo_[wh][st][se] = maxdrift;
0042   }
0043 }
0044 
0045 MuonPathConfirmator::~MuonPathConfirmator() {
0046   if (debug_)
0047     LogDebug("MuonPathConfirmator") << "MuonPathAnalyzer: destructor";
0048 }
0049 
0050 // ============================================================================
0051 // Main methods (initialise, run, finish)
0052 // ============================================================================
0053 
0054 void MuonPathConfirmator::run(edm::Event &iEvent,
0055                               const edm::EventSetup &iEventSetup,
0056                               std::vector<cmsdt::metaPrimitive> inMetaPrimitives,
0057                               edm::Handle<DTDigiCollection> dtdigis,
0058                               std::vector<cmsdt::metaPrimitive> &outMetaPrimitives) {
0059   if (debug_)
0060     LogDebug("MuonPathConfirmator") << "MuonPathConfirmator: run";
0061 
0062   // fit per SL (need to allow for multiple outputs for a single mpath)
0063   if (!inMetaPrimitives.empty()) {
0064     int dum_sl_rawid = inMetaPrimitives[0].rawId;
0065     DTSuperLayerId dumSlId(dum_sl_rawid);
0066     DTChamberId ChId(dumSlId.wheel(), dumSlId.station(), dumSlId.sector());
0067     max_drift_tdc = maxdriftinfo_[dumSlId.wheel() + 2][dumSlId.station() - 1][dumSlId.sector() - 1];
0068   }
0069 
0070   for (auto &mp : inMetaPrimitives) {
0071     analyze(mp, dtdigis, outMetaPrimitives);
0072   }
0073 }
0074 
0075 void MuonPathConfirmator::initialise(const edm::EventSetup &iEventSetup) {
0076   if (debug_)
0077     LogDebug("MuonPathConfirmator") << "MuonPathConfirmator::initialiase";
0078 }
0079 
0080 void MuonPathConfirmator::finish() {
0081   if (debug_)
0082     LogDebug("MuonPathConfirmator") << "MuonPathConfirmator: finish";
0083 };
0084 
0085 //------------------------------------------------------------------
0086 //--- Metodos privados
0087 //------------------------------------------------------------------
0088 
0089 void MuonPathConfirmator::analyze(cmsdt::metaPrimitive mp,
0090                                   edm::Handle<DTDigiCollection> dtdigis,
0091                                   std::vector<cmsdt::metaPrimitive> &outMetaPrimitives) {
0092   int dum_sl_rawid = mp.rawId;
0093   DTSuperLayerId dumSlId(dum_sl_rawid);
0094   DTChamberId ChId(dumSlId.wheel(), dumSlId.station(), dumSlId.sector());
0095   DTSuperLayerId sl1Id(ChId.rawId(), 1);
0096   DTSuperLayerId sl3Id(ChId.rawId(), 3);
0097 
0098   DTWireId wireIdSL1(sl1Id, 2, 1);
0099   DTWireId wireIdSL3(sl3Id, 2, 1);
0100   auto sl_shift_cm = shiftinfo_[wireIdSL1.rawId()] - shiftinfo_[wireIdSL3.rawId()];
0101   bool isSL1 = (mp.rawId == sl1Id.rawId());
0102   bool isSL3 = (mp.rawId == sl3Id.rawId());
0103   if (!isSL1 && !isSL3)
0104     outMetaPrimitives.emplace_back(mp);
0105   else {
0106     int best_tdc = -1;
0107     int next_tdc = -1;
0108     int best_wire = -1;
0109     int next_wire = -1;
0110     int best_layer = -1;
0111     int next_layer = -1;
0112     int best_lat = -1;
0113     int next_lat = -1;
0114     int lat = -1;
0115     int matched_digis = 0;
0116 
0117     int position_prec = ((int)(mp.x)) << PARTIALS_PRECISSION;
0118     int slope_prec = ((int)(mp.tanPhi)) << PARTIALS_PRECISSION;
0119 
0120     int slope_x_halfchamb = (((long int)slope_prec) * SEMICHAMBER_H) >> SEMICHAMBER_RES_SHR;
0121     int slope_x_3semicells = (slope_prec * 3) >> LYRANDAHALF_RES_SHR;
0122     int slope_x_1semicell = (slope_prec * 1) >> LYRANDAHALF_RES_SHR;
0123 
0124     for (const auto &dtLayerId_It : *dtdigis) {
0125       const DTLayerId dtLId = dtLayerId_It.first;
0126       // creating a new DTSuperLayerId object to compare with the required SL id
0127       const DTSuperLayerId dtSLId(dtLId.wheel(), dtLId.station(), dtLId.sector(), dtLId.superLayer());
0128       bool hitFromSL1 = (dtSLId.rawId() == sl1Id.rawId());
0129       bool hitFromSL3 = (dtSLId.rawId() == sl3Id.rawId());
0130       if (!(hitFromSL1 || hitFromSL3))  // checking hits are from one of the other SL of the same chamber
0131         continue;
0132       double minx = 10 * minx_match_2digis_ * ((double)max_drift_tdc / (double)CELL_SEMILENGTH);
0133       double min2x = 10 * minx_match_2digis_ * ((double)max_drift_tdc / (double)CELL_SEMILENGTH);
0134       if (isSL1 != hitFromSL1) {  // checking hits have the opposite SL than the TP
0135         for (auto digiIt = (dtLayerId_It.second).first; digiIt != (dtLayerId_It.second).second; ++digiIt) {
0136           if ((*digiIt).time() < mp.t0)
0137             continue;
0138           int wp_semicells = ((*digiIt).wire() - 1 - SL1_CELLS_OFFSET) * 2 + 1;
0139           int ly = dtLId.layer() - 1;
0140           if (ly % 2 == 1)
0141             wp_semicells -= 1;
0142           if (hitFromSL3)
0143             wp_semicells -= (int)round((sl_shift_cm * 10) / CELL_SEMILENGTH);
0144           double hit_position = wp_semicells * max_drift_tdc +
0145                                 ((*digiIt).time() - mp.t0) * (double)TIME_TO_TDC_COUNTS / (double)LHC_CLK_FREQ;
0146           double hit_position_left = wp_semicells * max_drift_tdc -
0147                                      ((*digiIt).time() - mp.t0) * (double)TIME_TO_TDC_COUNTS / (double)LHC_CLK_FREQ;
0148           // extrapolating position to the layer of the hit
0149           // mp.position is referred to the center between SLs, so one has to add half the distance between SLs
0150           // + half a cell height to get to the first wire + ly * cell height to reach the desired ly
0151           // 10 * VERT_PHI1_PHI3 / 2 + (CELL_HEIGHT / 2) + ly * CELL_HEIGHT = (10 * VERT_PHI1_PHI3 + (2 * ly + 1) * CELL_HEIGHT) / 2
0152 
0153           int position_in_layer = position_prec + (1 - 2 * (int)hitFromSL1) * slope_x_halfchamb;
0154           if (ly == 0)
0155             position_in_layer -= slope_x_3semicells;
0156           if (ly == 1)
0157             position_in_layer -= slope_x_1semicell;
0158           if (ly == 2)
0159             position_in_layer += slope_x_1semicell;
0160           if (ly == 3)
0161             position_in_layer += slope_x_3semicells;
0162           position_in_layer = position_in_layer >> PARTIALS_PRECISSION;
0163 
0164           if (std::abs(position_in_layer - hit_position_left) < std::abs(position_in_layer - hit_position)) {
0165             lat = 0;
0166             hit_position = hit_position_left;
0167           }
0168           if (std::abs(position_in_layer - hit_position) < minx) {
0169             // different layer than the stored in best, hit added, matched_digis++;. This approach in somewhat
0170             // buggy, as we could have stored as best LayerX -> LayerY -> LayerX, and this should
0171             // count only as 2 hits. However, as we confirm with at least 2 hits, having 2 or more
0172             // makes no difference
0173             if (dtLId.layer() != best_layer) {
0174               minx = std::abs(position_in_layer - hit_position);
0175               next_wire = best_wire;
0176               next_tdc = best_tdc;
0177               next_layer = best_layer;
0178               next_lat = best_lat;
0179               matched_digis++;
0180             }
0181             best_wire = (*digiIt).wire() - 1;
0182             best_tdc = (*digiIt).time();
0183             best_layer = dtLId.layer();
0184             best_lat = lat;
0185 
0186           } else if ((std::abs(position_in_layer - hit_position) >= minx) &&
0187                      (std::abs(position_in_layer - hit_position) < min2x)) {
0188             // same layer than the stored in best, no hit added
0189             if (dtLId.layer() == best_layer)
0190               continue;
0191             // different layer than the stored in next, hit added. This approach in somewhat
0192             // buggy, as we could have stored as next LayerX -> LayerY -> LayerX, and this should
0193             // count only as 2 hits. However, as we confirm with at least 2 hits, having 2 or more
0194             // makes no difference
0195             matched_digis++;
0196             // whether the layer is the same for this hit and the stored in next, we substitute
0197             // the one stored and modify the min distance
0198             min2x = std::abs(position_in_layer - hit_position);
0199             next_wire = (*digiIt).wire() - 1;
0200             next_tdc = (*digiIt).time();
0201             next_layer = dtLId.layer();
0202             next_lat = lat;
0203           }
0204         }
0205       }
0206     }
0207     int new_quality = mp.quality;
0208     std::vector<int> wi_c(4, -1), tdc_c(4, -1), lat_c(4, -1);
0209     if (matched_digis >= 2 and best_layer != -1 and next_layer != -1) {  // actually confirm
0210       new_quality = CHIGHQ;
0211       if (mp.quality == LOWQ)
0212         new_quality = CLOWQ;
0213 
0214       wi_c[next_layer - 1] = next_wire;
0215       tdc_c[next_layer - 1] = next_tdc;
0216       lat_c[next_layer - 1] = next_lat;
0217 
0218       wi_c[best_layer - 1] = best_wire;
0219       tdc_c[best_layer - 1] = best_tdc;
0220       lat_c[best_layer - 1] = best_lat;
0221     }
0222     if (isSL1) {
0223       outMetaPrimitives.emplace_back(metaPrimitive(
0224           {mp.rawId, mp.t0,   mp.x,     mp.tanPhi, mp.phi,   mp.phiB, mp.phi_cmssw, mp.phiB_cmssw, mp.chi2, new_quality,
0225            mp.wi1,   mp.tdc1, mp.lat1,  mp.wi2,    mp.tdc2,  mp.lat2, mp.wi3,       mp.tdc3,       mp.lat3, mp.wi4,
0226            mp.tdc4,  mp.lat4, wi_c[0],  tdc_c[0],  lat_c[0], wi_c[1], tdc_c[1],     lat_c[1],      wi_c[2], tdc_c[2],
0227            lat_c[2], wi_c[3], tdc_c[3], lat_c[3],  -1}));
0228     } else {
0229       outMetaPrimitives.emplace_back(
0230           metaPrimitive({mp.rawId,      mp.t0,    mp.x,        mp.tanPhi, mp.phi,   mp.phiB,  mp.phi_cmssw,
0231                          mp.phiB_cmssw, mp.chi2,  new_quality, wi_c[0],   tdc_c[0], lat_c[0], wi_c[1],
0232                          tdc_c[1],      lat_c[1], wi_c[2],     tdc_c[2],  lat_c[2], wi_c[3],  tdc_c[3],
0233                          lat_c[3],      mp.wi5,   mp.tdc5,     mp.lat5,   mp.wi6,   mp.tdc6,  mp.lat6,
0234                          mp.wi7,        mp.tdc7,  mp.lat7,     mp.wi8,    mp.tdc8,  mp.lat8,  -1}));
0235     }
0236   }  //SL2
0237 }