DTTriggerPhase2/src/RPCIntegrator.cc

0001 #include "FWCore/Utilities/interface/ESGetToken.h"
0002 #include "FWCore/Utilities/interface/ESInputTag.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "L1Trigger/DTTriggerPhase2/interface/RPCIntegrator.h"
0006
0007 #include <Geometry/Records/interface/MuonGeometryRecord.h>
0008
0009 #include <cmath>
0010
0011 using namespace cmsdt;
0012
0013 RPCIntegrator::RPCIntegrator(const edm::ParameterSet& pset, edm::ConsumesCollector& iC)
0014     : m_debug_(pset.getUntrackedParameter<bool>("debug")),
0015       m_max_quality_to_overwrite_t0_(pset.getParameter<int>("max_quality_to_overwrite_t0")),
0016       m_bx_window_(pset.getParameter<int>("bx_window")),
0017       m_phi_window_(pset.getParameter<double>("phi_window")),
0018       m_storeAllRPCHits_(pset.getParameter<bool>("storeAllRPCHits")) {
0019   if (m_debug_)
0020     LogDebug("RPCIntegrator") << "RPCIntegrator constructor";
0021
0022   rpcGeomH_ = iC.esConsumes<RPCGeometry, MuonGeometryRecord>();
0023   dtGeomH_ = iC.esConsumes<DTGeometry, MuonGeometryRecord>();
0024 }
0025
0026 RPCIntegrator::~RPCIntegrator() {
0027   if (m_debug_)
0028     LogDebug("RPCIntegrator") << "RPCIntegrator destructor";
0029 }
0030
0031 void RPCIntegrator::initialise(const edm::EventSetup& iEventSetup, double shift_back_fromDT) {
0032   if (m_debug_)
0033     LogDebug("RPCIntegrator") << "RPCIntegrator initialisation";
0034
0035   if (m_debug_)
0036     LogDebug("RPCIntegrator") << "Getting RPC geometry";
0037
0038   if (auto handle = iEventSetup.getHandle(dtGeomH_)) {
0039     dtGeo_ = handle.product();
0040   }
0041
0042   if (auto handle = iEventSetup.getHandle(rpcGeomH_)) {
0043     rpcGeo_ = handle.product();
0044   }
0045
0046   shift_back_ = shift_back_fromDT;
0047 }
0048
0049 void RPCIntegrator::finish() {}
0050
0051 void RPCIntegrator::prepareMetaPrimitives(edm::Handle<RPCRecHitCollection> rpcRecHits) {
0052   RPCMetaprimitives_.clear();
0053   rpcRecHits_translated_.clear();
0054   for (const auto& rpcIt : *rpcRecHits) {
0055     RPCDetId rpcDetId = (RPCDetId)(rpcIt).rpcId();
0056     GlobalPoint global_position = RPCGlobalPosition(rpcDetId, rpcIt);
0057     int rpc_region = rpcDetId.region();
0058     if (rpc_region != 0)
0059       continue;  // Region = 0 Barrel
0060
0061     // set everyone to rpc single hit (3) not matched to DT flag for now
0062     // change last two elements if dt bx centered at zero again
0063     RPCMetaprimitives_.emplace_back(
0064         rpcDetId, &rpcIt, global_position, RPC_HIT, rpcIt.BunchX() + BX_SHIFT, rpcIt.time() + BX_SHIFT * LHC_CLK_FREQ);
0065   }
0066 }
0067 void RPCIntegrator::matchWithDTAndUseRPCTime(std::vector<metaPrimitive>& dt_metaprimitives) {
0068   for (auto dt_metaprimitive = dt_metaprimitives.begin(); dt_metaprimitive != dt_metaprimitives.end();
0069        dt_metaprimitive++) {
0070     RPCMetaprimitive* bestMatch_rpcRecHit = matchDTwithRPC(&*dt_metaprimitive);
0071     if (bestMatch_rpcRecHit) {
0072       (*dt_metaprimitive).rpcFlag = RPC_CONFIRM;
0073       if ((*dt_metaprimitive).quality < m_max_quality_to_overwrite_t0_) {
0074         (*dt_metaprimitive).t0 = bestMatch_rpcRecHit->rpc_t0 + 25 * shift_back_;
0075         (*dt_metaprimitive).rpcFlag = RPC_TIME;
0076       }
0077     }
0078   }
0079 }
0080
0081 void RPCIntegrator::makeRPCOnlySegments() {
0082   std::vector<L1Phase2MuDTPhDigi> rpc_only_segments;
0083   for (auto& rpc_mp_it_layer1 : RPCMetaprimitives_) {
0084     RPCDetId rpc_id_l1 = rpc_mp_it_layer1.rpc_id;
0085     const RPCRecHit* rpc_cluster_l1 = rpc_mp_it_layer1.rpc_cluster;
0086     GlobalPoint rpc_gp_l1 = rpc_mp_it_layer1.global_position;
0087     if (rpc_id_l1.station() > 2 || rpc_id_l1.layer() != 1 ||
0088         (rpc_mp_it_layer1.rpcFlag == RPC_ASSOCIATE && !m_storeAllRPCHits_))
0089       continue;
0090     // only one RPC layer in station three and four &&
0091     // avoid duplicating pairs &&
0092     // avoid building RPC only segment if DT segment was already there
0093     int min_dPhi = std::numeric_limits<int>::max();
0094     RPCMetaprimitive* bestMatch_rpc_mp_layer2 = nullptr;
0095     for (auto& rpc_mp_it_layer2 : RPCMetaprimitives_) {
0096       RPCDetId rpc_id_l2 = rpc_mp_it_layer2.rpc_id;
0097       const RPCRecHit* rpc_cluster_l2 = rpc_mp_it_layer2.rpc_cluster;
0098       GlobalPoint rpc_gp_l2 = rpc_mp_it_layer2.global_position;
0099       if (rpc_id_l2.station() == rpc_id_l1.station() && rpc_id_l2.ring() == rpc_id_l1.ring() &&
0100           rpc_id_l2.layer() != rpc_id_l1.layer()  // ensure to have layer 1 --> layer 2
0101           && rpc_id_l2.sector() == rpc_id_l1.sector() && rpc_cluster_l2->BunchX() == rpc_cluster_l1->BunchX() &&
0102           (rpc_mp_it_layer2.rpcFlag != RPC_ASSOCIATE || m_storeAllRPCHits_)) {
0103         // avoid building RPC only segment with a hit already matched to DT,
0104         // except if one aske to store all RPC info
0105         float tmp_dPhi = rpc_gp_l1.phi() - rpc_gp_l2.phi();
0106         if (std::abs(tmp_dPhi) < std::abs(min_dPhi)) {
0107           min_dPhi = tmp_dPhi;
0108           bestMatch_rpc_mp_layer2 = &rpc_mp_it_layer2;
0109         }
0110       }
0111     }
0112     if (bestMatch_rpc_mp_layer2) {
0113       rpc_mp_it_layer1.rpcFlag = 6;
0114       // need a new flag (will be removed later) to differentiate
0115       // between "has been matched to DT" and "Has been used in an RPC only segment"
0116       bestMatch_rpc_mp_layer2->rpcFlag = 6;
0117       double phiB = phiBending(&rpc_mp_it_layer1, &*bestMatch_rpc_mp_layer2);
0118       // Arbitrarily choose the phi from layer 1
0119       double global_phi = rpc_mp_it_layer1.global_position.phi();
0120       double t0 = (rpc_mp_it_layer1.rpc_t0 + bestMatch_rpc_mp_layer2->rpc_t0) / 2;
0121       // RPC only segment have rpcFlag==2
0122       L1Phase2MuDTPhDigi rpc_only_segment =
0123           createL1Phase2MuDTPhDigi(rpc_id_l1, rpc_mp_it_layer1.rpc_bx, t0, global_phi, phiB, 2);
0124       rpc_only_segments.push_back(rpc_only_segment);
0125     }
0126   }
0127   rpcRecHits_translated_.insert(rpcRecHits_translated_.end(), rpc_only_segments.begin(), rpc_only_segments.end());
0128 }
0129
0130 void RPCIntegrator::storeRPCSingleHits() {
0131   for (auto rpc_mp_it = RPCMetaprimitives_.begin(); rpc_mp_it != RPCMetaprimitives_.end(); rpc_mp_it++) {
0132     RPCDetId rpcDetId = rpc_mp_it->rpc_id;
0133     if (rpc_mp_it->rpcFlag == 6)
0134       rpc_mp_it->rpcFlag = RPC_ASSOCIATE;
0135     L1Phase2MuDTPhDigi rpc_out = createL1Phase2MuDTPhDigi(
0136         rpcDetId, rpc_mp_it->rpc_bx, rpc_mp_it->rpc_t0, rpc_mp_it->global_position.phi(), -10000, rpc_mp_it->rpcFlag);
0137     rpcRecHits_translated_.push_back(rpc_out);
0138   }
0139 }
0140
0141 void RPCIntegrator::removeRPCHitsUsed() {
0142   if (m_debug_)
0143     LogDebug("RPCIntegrator") << "RPCIntegrator removeRPCHitsUsed method";
0144   if (!m_storeAllRPCHits_) {
0145     // Remove RPC hit attached to a DT or RPC segment if required by user
0146     // (avoid having two TP's corresponding to the same physical hit)
0147     auto rpcRecHit_translated_ = rpcRecHits_translated_.begin();
0148     while (rpcRecHit_translated_ != rpcRecHits_translated_.end()) {
0149       if (rpcRecHit_translated_->rpcFlag() == RPC_ASSOCIATE || rpcRecHit_translated_->rpcFlag() == 6) {
0150         rpcRecHit_translated_ = rpcRecHits_translated_.erase(rpcRecHit_translated_);
0151       } else {
0152         ++rpcRecHit_translated_;
0153       }
0154     }
0155   }
0156 }
0157
0158 RPCMetaprimitive* RPCIntegrator::matchDTwithRPC(metaPrimitive* dt_metaprimitive) {
0159   // metaprimitive dtChId is still in convention with [1 - 12]
0160   // because at this stage the BX of metaprimitive is not yet computed...
0161   // will also have to subtract 20*25 ns because of the recent change
0162   int dt_bx = (int)round(dt_metaprimitive->t0 / 25.) - shift_back_;
0163   DTChamberId dt_chId = DTChamberId(dt_metaprimitive->rawId);
0164   int dt_sector = dt_chId.sector();
0165   if (dt_sector == 13)
0166     dt_sector = 4;
0167   if (dt_sector == 14)
0168     dt_sector = 10;
0169   RPCMetaprimitive* bestMatch_rpcRecHit = nullptr;
0170   float min_dPhi = std::numeric_limits<float>::max();
0171   for (auto rpc_mp_it = RPCMetaprimitives_.begin(); rpc_mp_it != RPCMetaprimitives_.end(); rpc_mp_it++) {
0172     RPCDetId rpc_det_id = rpc_mp_it->rpc_id;
0173     if (rpc_det_id.ring() == dt_chId.wheel()  // ring() in barrel RPC corresponds to the wheel
0174         && rpc_det_id.station() == dt_chId.station() && rpc_det_id.sector() == dt_sector &&
0175         std::abs(rpc_mp_it->rpc_bx - dt_bx) <= m_bx_window_) {
0176       // Select the RPC hit closest in phi to the DT meta primitive
0177
0178       // just a trick to apply the phi window cut on what could be accessed to fine tune it
0179       int delta_phi =
0180           (int)round((phi_DT_MP_conv(rpc_mp_it->global_position.phi(), rpc_det_id.sector()) - dt_metaprimitive->phi) *
0181                      cmsdt::PHIBRES_CONV);
0182       if (std::abs(delta_phi) < min_dPhi && std::abs(delta_phi) < m_phi_window_) {
0183         min_dPhi = std::abs(delta_phi);
0184         bestMatch_rpcRecHit = &*rpc_mp_it;
0185       }
0186     }
0187   }
0188   if (bestMatch_rpcRecHit) {
0189     bestMatch_rpcRecHit->rpcFlag = RPC_ASSOCIATE;
0190   }
0191   return bestMatch_rpcRecHit;
0192 }
0193
0194 L1Phase2MuDTPhDigi RPCIntegrator::createL1Phase2MuDTPhDigi(
0195     RPCDetId rpcDetId, int rpc_bx, double rpc_time, double rpc_global_phi, double phiB, int rpc_flag) {
0196   if (m_debug_)
0197     LogDebug("RPCIntegrator") << "Creating DT TP out of RPC recHits";
0198   int rpc_wheel = rpcDetId.ring();             // In barrel, wheel is accessed via ring() method ([-2,+2])
0199   int trigger_sector = rpcDetId.sector() - 1;  // DT Trigger sector:[0,11] while RPC sector:[1,12]
0200   int rpc_station = rpcDetId.station();
0201   int rpc_layer = rpcDetId.layer();
0202   int rpc_trigger_phi = phiInDTTPFormat(rpc_global_phi, rpcDetId.sector());
0203   int rpc_trigger_phiB = (phiB == -10000) ? phiB : (int)round(phiB * cmsdt::PHIBRES_CONV);
0204   int rpc_quality = -1;  // dummy for rpc
0205   int rpc_index = 0;     // dummy for rpc
0206   return L1Phase2MuDTPhDigi(rpc_bx,
0207                             rpc_wheel,
0208                             trigger_sector,
0209                             rpc_station,
0210                             rpc_layer,  //this would be the layer in the new dataformat
0211                             rpc_trigger_phi,
0212                             rpc_trigger_phiB,
0213                             rpc_quality,
0214                             rpc_index,
0215                             rpc_time,
0216                             -1,  // no chi2 for RPC
0217                             rpc_flag);
0218 }
0219
0220 double RPCIntegrator::phiBending(RPCMetaprimitive* rpc_hit_1, RPCMetaprimitive* rpc_hit_2) {
0221   DTChamberId DT_chamber(rpc_hit_1->rpc_id.ring(), rpc_hit_1->rpc_id.station(), rpc_hit_1->rpc_id.sector());
0222   LocalPoint lp_rpc_hit_1_dtconv = dtGeo_->chamber(DT_chamber)->toLocal(rpc_hit_1->global_position);
0223   LocalPoint lp_rpc_hit_2_dtconv = dtGeo_->chamber(DT_chamber)->toLocal(rpc_hit_2->global_position);
0224   double slope = (lp_rpc_hit_1_dtconv.x() - lp_rpc_hit_2_dtconv.x()) / distance_between_two_rpc_layers_;
0225   double average_x = (lp_rpc_hit_1_dtconv.x() + lp_rpc_hit_2_dtconv.x()) / 2;
0226   GlobalPoint seg_middle_global =
0227       dtGeo_->chamber(DT_chamber)->toGlobal(LocalPoint(average_x, 0., 0.));  // for station 1 and 2, z = 0
0228   double seg_phi = phi_DT_MP_conv(seg_middle_global.phi(), rpc_hit_1->rpc_id.sector());
0229   double psi = atan(slope);
0230   double phiB = hasPosRF_rpc(rpc_hit_1->rpc_id.ring(), rpc_hit_1->rpc_id.sector()) ? psi - seg_phi : -psi - seg_phi;
0231   return phiB;
0232 }
0233
0234 int RPCIntegrator::phiInDTTPFormat(double rpc_global_phi, int rpcSector) {
0235   double rpc_localDT_phi;
0236   rpc_localDT_phi = phi_DT_MP_conv(rpc_global_phi, rpcSector) * cmsdt::PHIBRES_CONV;
0237   return (int)round(rpc_localDT_phi);
0238 }
0239
0240 double RPCIntegrator::phi_DT_MP_conv(double rpc_global_phi, int rpcSector) {
0241   // Adaptation of https://github.com/cms-sw/cmssw/blob/master/L1Trigger/L1TTwinMux/src/RPCtoDTTranslator.cc#L349
0242
0243   if (rpcSector == 1)
0244     return rpc_global_phi;
0245   else {
0246     float conversion = 1 / 6.;
0247     if (rpc_global_phi >= 0)
0248       return rpc_global_phi - (rpcSector - 1) * M_PI * conversion;
0249     else
0250       return rpc_global_phi + (13 - rpcSector) * M_PI * conversion;
0251   }
0252 }
0253
0254 GlobalPoint RPCIntegrator::RPCGlobalPosition(RPCDetId rpcId, const RPCRecHit& rpcIt) const {
0255   RPCDetId rpcid = RPCDetId(rpcId);
0256   const LocalPoint& rpc_lp = rpcIt.localPosition();
0257   const GlobalPoint& rpc_gp = rpcGeo_->idToDet(rpcid)->surface().toGlobal(rpc_lp);
0258
0259   return rpc_gp;
0260 }
0261
0262 bool RPCIntegrator::hasPosRF_rpc(int wh, int sec) const { return (wh > 0 || (wh == 0 && sec % 4 > 1)); }