#include <AngleCalculation.h>

Public Member Functions
void	calculate_angles (EMTFTrack &track) const

void	calculate_bx (EMTFTrack &track) const

void	configure (int verbose, int endcap, int sector, int bx, int bxWindow, int thetaWindow, int thetaWindowRPC, bool bugME11Dupes)

void	erase_tracks (EMTFTrackCollection &tracks) const

void	process (emtf::zone_array< EMTFTrackCollection > &zone_tracks) const

Private Attributes
bool	bugME11Dupes_

int	bx_

int	bxWindow_

int	endcap_

int	sector_

int	thetaWindow_

int	thetaWindowRPC_

int	verbose_

Detailed Description

Definition at line 7 of file AngleCalculation.h.

Member Function Documentation

void AngleCalculation::calculate_angles ( EMTFTrack & track ) const

Definition at line 82 of file AngleCalculation.cc.

Referenced by process().

                                                               {
   // Group track hits by station
   std::array<EMTFHitCollection, emtf::NUM_STATIONS> st_conv_hits;
 
   for (int istation = 0; istation < emtf::NUM_STATIONS; ++istation) {
     for (const auto& conv_hit : track.Hits()) {
       if ((conv_hit.Station() - 1) == istation) {
         st_conv_hits.at(istation).push_back(conv_hit);
       }
     }
 
     if (bugME11Dupes_) {
       if (not (st_conv_hits.at(istation).size() <= 4)) // ambiguity in theta is max 4
     { edm::LogError("L1T") << "st_conv_hits.at(istation).size() = " << st_conv_hits.at(istation).size(); return; }
     } else {
       if (not (st_conv_hits.at(istation).size() <= 2))  // ambiguity in theta is max 2
     { edm::LogError("L1T") << "st_conv_hits.at(istation).size() = " << st_conv_hits.at(istation).size(); return; }
     }
   }
   if (not (st_conv_hits.size() == emtf::NUM_STATIONS))
     { edm::LogError("L1T") << "st_conv_hits.size() = " << st_conv_hits.size() << ", emtf::NUM_STATIONS = " << emtf::NUM_STATIONS; return; }
 
 
   // Best theta deltas and phi deltas
   // from 0 to 5: dtheta12, dtheta13, dtheta14, dtheta23, dtheta24, dtheta34
   std::array<int,  emtf::NUM_STATION_PAIRS> best_dtheta_arr;  // Best of up to 4 dTheta values per pair of stations (with duplicate thetas)
   std::array<int,  emtf::NUM_STATION_PAIRS> best_dtheta_sign_arr;
   std::array<int,  emtf::NUM_STATION_PAIRS> best_dphi_arr;   // Not really "best" - there is only one dPhi value per pair of stations
   std::array<int,  emtf::NUM_STATION_PAIRS> best_dphi_sign_arr;
 
   // Best angles
   // from 0 to 5: ME2,      ME3,      ME4,      ME2,      ME2,      ME3
   //              dtheta12, dtheta13, dtheta14, dtheta23, dtheta24, dtheta34
   std::array<int,  emtf::NUM_STATION_PAIRS> best_theta_arr;
   std::array<int,  emtf::NUM_STATION_PAIRS> best_phi_arr;
 
   // Keep track of which pair is valid
   std::array<bool, emtf::NUM_STATION_PAIRS> best_dtheta_valid_arr;
   std::array<bool, emtf::NUM_STATION_PAIRS> best_has_rpc_arr; // Not used - should remove (AWB 21.05.17)
 
   // Initialize
   best_dtheta_arr      .fill(invalid_dtheta);
   best_dtheta_sign_arr .fill(1);
   best_dphi_arr        .fill(invalid_dphi);
   best_dphi_sign_arr   .fill(1);
   best_phi_arr         .fill(0);
   best_theta_arr       .fill(0);
   best_dtheta_valid_arr.fill(false);
   best_has_rpc_arr     .fill(false);
 
   auto abs_diff = [](int a, int b) { return std::abs(a-b); };
 
   // Calculate angles
   int ipair = 0;
 
   for (int ist1 = 0; ist1+1 < emtf::NUM_STATIONS; ++ist1) {  // station A
     for (int ist2 = ist1+1; ist2 < emtf::NUM_STATIONS; ++ist2) {  // station B
       const EMTFHitCollection& conv_hitsA = st_conv_hits.at(ist1);
       const EMTFHitCollection& conv_hitsB = st_conv_hits.at(ist2);
 
       // More than 1 hit per station when hit has ambigous theta
       for (const auto& conv_hitA : conv_hitsA) {
         for (const auto& conv_hitB : conv_hitsB) {
           // Has RPC?
           //bool has_rpc = (conv_hitA.Subsystem() == TriggerPrimitive::kRPC || conv_hitB.Subsystem() == TriggerPrimitive::kRPC);
 
           // Calculate theta deltas
           int thA = conv_hitA.Theta_fp();
           int thB = conv_hitB.Theta_fp();
           int dth = abs_diff(thA, thB);
           int dth_sign = (thA <= thB);  // sign
           if (not(thA != 0 && thB != 0))
         { edm::LogError("L1T") << "thA = " << thA << ", thB = " << thB; return; }
           if (not(dth < invalid_dtheta))
         { edm::LogError("L1T") << "dth = " << dth << ", invalid_dtheta = " << invalid_dtheta; return; }
 
           if (best_dtheta_arr.at(ipair) >= dth) {
             best_dtheta_arr.at(ipair) = dth;
             best_dtheta_sign_arr.at(ipair) = dth_sign;
             best_dtheta_valid_arr.at(ipair) = true;
             // best_has_rpc_arr.at(ipair) = has_rpc;  // FW doesn't check whether a segment is CSC or RPC, should remove - AWB 21.05.17
 
             // first 3 pairs, use station B
             // last 3 pairs, use station A
             best_theta_arr.at(ipair) = (ipair < 3) ? thB : thA;
           }
 
           // Calculate phi deltas
           int phA = conv_hitA.Phi_fp();
           int phB = conv_hitB.Phi_fp();
           int dph = abs_diff(phA, phB);
           int dph_sign = (phA <= phB);
 
           if (best_dphi_arr.at(ipair) >= dph) {
             best_dphi_arr.at(ipair) = dph;
             best_dphi_sign_arr.at(ipair) = dph_sign;
 
             // first 3 pairs, use station B
             // last 3 pairs, use station A
             best_phi_arr.at(ipair) = (ipair < 3) ? phB : phA;
           }
         }  // end loop over conv_hits in station B
       }  // end loop over conv_hits in station A
 
       ++ipair;
     }  // end loop over station B
   }  // end loop over station A
   if (not(ipair == emtf::NUM_STATION_PAIRS))
     { edm::LogError("L1T") << "ipair = " << ipair << ", emtf::NUM_STATION_PAIRS = " << emtf::NUM_STATION_PAIRS; return; }
 
 
   // Apply cuts on dtheta
 
   // There is a possible bug in FW. After a dtheta pair fails the theta window
   // cut, the valid flag of the pair is not updated. Later on, theta from
   // this pair is used to assign the precise theta of the track.
   std::array<bool, emtf::NUM_STATION_PAIRS> best_dtheta_valid_arr_1;
 
   for (int ipair = 0; ipair < emtf::NUM_STATION_PAIRS; ++ipair) {
     if (best_has_rpc_arr.at(ipair))
       best_dtheta_valid_arr_1.at(ipair) = best_dtheta_valid_arr.at(ipair) && (best_dtheta_arr.at(ipair) <= thetaWindowRPC_);
     else
       best_dtheta_valid_arr_1.at(ipair) = best_dtheta_valid_arr.at(ipair) && (best_dtheta_arr.at(ipair) <= thetaWindow_);
   }
 
   // Find valid segments
   // vmask contains valid station mask = {ME4,ME3,ME2,ME1}. "0b" prefix for binary.
   int vmask1 = 0, vmask2 = 0, vmask3 = 0;
 
   if (best_dtheta_valid_arr_1.at(0)) {
     vmask1 |= 0b0011;  // 12
   }
   if (best_dtheta_valid_arr_1.at(1)) {
     vmask1 |= 0b0101;  // 13
   }
   if (best_dtheta_valid_arr_1.at(2)) {
     vmask1 |= 0b1001;  // 14
   }
   if (best_dtheta_valid_arr_1.at(3)) {
     vmask2 |= 0b0110;  // 23
   }
   if (best_dtheta_valid_arr_1.at(4)) {
     vmask2 |= 0b1010;  // 24
   }
   if (best_dtheta_valid_arr_1.at(5)) {
     vmask3 |= 0b1100;  // 34
   }
 
   // merge station masks only if they share bits
   // Station 1 hits pass if any dTheta1X values pass
   // Station 2 hits pass if any dTheta2X values pass, *EXCEPT* the following cases:
   //           Only {13, 24} pass, only {13, 24, 34} pass,
   //           Only {14, 23} pass, only {14, 23, 34} pass.
   // Station 3 hits pass if any dTheta3X values pass, *EXCEPT* the following cases:
   //           Only {12, 34} pass, only {14, 23} pass.
   // Station 4 hits pass if any dTheta4X values pass, *EXCEPT* the following cases:
   //           Only {12, 34} pass, only {13, 24} pass.
   int vstat = vmask1;  // valid stations based on th coordinates
   if ((vstat & vmask2) != 0 || vstat == 0)
     vstat |= vmask2;
   if ((vstat & vmask3) != 0 || vstat == 0)
     vstat |= vmask3;
 
   // remove valid flag for station if hit does not pass the dTheta mask
   for (int istation = 0; istation < emtf::NUM_STATIONS; ++istation) {
     if ((vstat & (1 << istation)) == 0) {  // station bit not set
       st_conv_hits.at(istation).clear();
     }
   }
 
   // assign precise phi and theta for the track
   int phi_fp    = 0;
   int theta_fp  = 0;
   int best_pair = -1;
 
   if ((vstat & (1<<1)) != 0) {            // ME2 present
     if (!best_has_rpc_arr.at(0) && best_dtheta_valid_arr.at(0))      // 12, no RPC
       best_pair = 0;
     else if (!best_has_rpc_arr.at(3) && best_dtheta_valid_arr.at(3)) // 23, no RPC
       best_pair = 3;
     else if (!best_has_rpc_arr.at(4) && best_dtheta_valid_arr.at(4)) // 24, no RPC
       best_pair = 4;
     else if (best_dtheta_valid_arr.at(0)) // 12, has RPC
       best_pair = 0;
     else if (best_dtheta_valid_arr.at(3)) // 23, has RPC
       best_pair = 3;
     else if (best_dtheta_valid_arr.at(4)) // 24, has RPC
       best_pair = 4;
   } else if ((vstat & (1<<2)) != 0) {     // ME3 present
     if (!best_has_rpc_arr.at(1) && best_dtheta_valid_arr.at(1))      // 13, no RPC
       best_pair = 1;
     else if (!best_has_rpc_arr.at(5) && best_dtheta_valid_arr.at(5)) // 34, no RPC
       best_pair = 5;
     else if (best_dtheta_valid_arr.at(1)) // 13, has RPC
       best_pair = 1;
     else if (best_dtheta_valid_arr.at(5)) // 34, has RPC
       best_pair = 5;
   } else if ((vstat & (1<<3)) != 0) {     // ME4 present
     if (best_dtheta_valid_arr.at(2))      // 14
       best_pair = 2;
   }
 
   if (best_pair != -1) {
     phi_fp   = best_phi_arr.at(best_pair);
     theta_fp = best_theta_arr.at(best_pair);
     if (not(theta_fp != 0))
       { edm::LogError("L1T") << "theta_fp = " << theta_fp; return; }
 
     // In firmware, the track is associated to LCTs by the segment number, which
     // identifies the best strip, but does not resolve the ambiguity in theta.
     // In emulator, this additional logic also resolves the ambiguity in theta.
     struct {
       typedef EMTFHit value_type;
       bool operator()(const value_type& lhs, const value_type& rhs) const {
         return std::abs(lhs.Theta_fp()-theta) < std::abs(rhs.Theta_fp()-theta);
       }
       int theta;
     } less_dtheta_cmp;
     less_dtheta_cmp.theta = theta_fp;  // capture
 
     for (int istation = 0; istation < emtf::NUM_STATIONS; ++istation) {
       std::stable_sort(st_conv_hits.at(istation).begin(), st_conv_hits.at(istation).end(), less_dtheta_cmp);
       if (st_conv_hits.at(istation).size() > 1)
         st_conv_hits.at(istation).resize(1);  // just the minimum in dtheta
     }
   }
 
   // update rank taking into account available stations after theta deltas
   // keep straightness as it was
   int old_rank = (track.Rank() << 1);  // output rank is one bit longer than input rank, to accomodate ME4 separately
   int rank = (
       (((old_rank >> 6) & 1) << 6) |  // straightness
       (((old_rank >> 4) & 1) << 4) |  // straightness
       (((old_rank >> 2) & 1) << 2) |  // straightness
       (((vstat >> 0) & 1) << 5) |  // ME1
       (((vstat >> 1) & 1) << 3) |  // ME2
       (((vstat >> 2) & 1) << 1) |  // ME3
       (((vstat >> 3) & 1) << 0)    // ME4
   );
 
   int mode = (
       (((vstat >> 0) & 1) << 3) |  // ME1
       (((vstat >> 1) & 1) << 2) |  // ME2
       (((vstat >> 2) & 1) << 1) |  // ME3
       (((vstat >> 3) & 1) << 0)    // ME4
   );
 
   int mode_inv = vstat;
 
   // if less than 2 segments, kill rank
   if (vstat == 0b0001 || vstat == 0b0010 || vstat == 0b0100 || vstat == 0b1000 || vstat == 0)
     rank = 0;
 
   // From RecoMuon/DetLayers/src/MuonCSCDetLayerGeometryBuilder.cc
   auto isFront = [](int station, int ring, int chamber, int subsystem) {
 
     // // RPCs are behind the CSCs in stations 1, 3, and 4; in front in 2
     // if (subsystem == TriggerPrimitive::kRPC)
     //   return (station == 2);
 
     // In EMTF firmware, RPC hits are treated as if they came from the corresponding 
     // CSC chamber, so the FR bit assignment is the same as for CSCs - AWB 06.06.17
 
     // GEMs are in front of the CSCs
     if (subsystem == TriggerPrimitive::kGEM)
       return true;
 
     bool result = false;
     bool isOverlapping = !(station == 1 && ring == 3);
     // not overlapping means back
     if(isOverlapping)
     {
       bool isEven = (chamber % 2 == 0);
       // odd chambers are bolted to the iron, which faces
       // forward in 1&2, backward in 3&4, so...
       result = (station < 3) ? isEven : !isEven;
     }
     return result;
   };
 
   // Fill ptlut_data
   EMTFPtLUT ptlut_data = {};
   for (int i = 0; i < emtf::NUM_STATION_PAIRS; ++i) {
     ptlut_data.delta_ph[i] = best_dphi_arr.at(i);
     ptlut_data.sign_ph[i]  = best_dphi_sign_arr.at(i);
     ptlut_data.delta_th[i] = best_dtheta_arr.at(i);
     ptlut_data.sign_th[i]  = best_dtheta_sign_arr.at(i);
   }
 
   for (int i = 0; i < emtf::NUM_STATIONS; ++i) {
     const auto& v = st_conv_hits.at(i);
     ptlut_data.cpattern[i] = v.empty() ? 0 : v.front().Pattern();  // Automatically set to 0 for RPCs
     ptlut_data.fr[i]       = v.empty() ? 0 : isFront(v.front().Station(), v.front().Ring(), v.front().Chamber(), v.front().Subsystem());
     if (i == 0) 
       ptlut_data.st1_ring2 = v.empty() ? 0 : (v.front().Station() == 1 && (v.front().Ring() == 2 || v.front().Ring() == 3));
   }
 
   for (int i = 0; i < emtf::NUM_STATIONS+1; ++i) {  // 'bt' arrays use 5-station convention
     ptlut_data.bt_vi[i] = 0;
     ptlut_data.bt_hi[i] = 0;
     ptlut_data.bt_ci[i] = 0;
     ptlut_data.bt_si[i] = 0;
   }
 
   for (int i = 0; i < emtf::NUM_STATIONS; ++i) {
     const auto& v = st_conv_hits.at(i);
     if (!v.empty()) {
       int bt_station = v.front().BT_station();
       if (not(0 <= bt_station && bt_station <= 4))
     { edm::LogError("L1T") << "bt_station = " << bt_station; return; }
 
       int bt_segment = v.front().BT_segment();
       ptlut_data.bt_vi[bt_station] = 1;
       ptlut_data.bt_hi[bt_station] = (bt_segment >> 5) & 0x3;
       ptlut_data.bt_ci[bt_station] = (bt_segment >> 1) & 0xf;
       ptlut_data.bt_si[bt_station] = (bt_segment >> 0) & 0x1;
     }
   }
 
   // ___________________________________________________________________________
   // Output
 
   track.set_rank     ( rank );
   track.set_mode     ( mode );
   track.set_mode_inv ( mode_inv );
   track.set_phi_fp   ( phi_fp );
   track.set_theta_fp ( theta_fp );
   track.set_PtLUT    ( ptlut_data );
 
   track.set_phi_loc  ( emtf::calc_phi_loc_deg(phi_fp) );
   track.set_phi_glob ( emtf::calc_phi_glob_deg(track.Phi_loc(), track.Sector()) );
   track.set_theta    ( emtf::calc_theta_deg_from_int(theta_fp) );
   track.set_eta      ( emtf::calc_eta_from_theta_deg(track.Theta(), track.Endcap()) );
 
   // Only keep the best segments
   track.clear_Hits();
 
   EMTFHitCollection tmp_hits = track.Hits();
   flatten_container(st_conv_hits, tmp_hits);
   track.set_Hits( tmp_hits );
 
 }

void AngleCalculation::calculate_bx ( EMTFTrack & track ) const

Definition at line 425 of file AngleCalculation.cc.

References bx_, bxWindow_, counter, l1t::EMTFTrack::Hits(), mps_fire::i, l1t::EMTFTrack::set_first_bx(), and l1t::EMTFTrack::set_second_bx().

Referenced by process().

                                                           {
   const int delayBX = bxWindow_ - 1;
   if (not(delayBX >= 0))
     { edm::LogError("L1T") << "delayBX = " << delayBX; return; }
   std::vector<int> counter(delayBX+1, 0);
 
   for (const auto& conv_hit : track.Hits()) {
     for (int i = delayBX; i >= 0; i--) {
       if (conv_hit.BX() <= bx_ - i)
         counter.at(i) += 1;  // Count stubs delayed by i BX or more
     }
   }
 
   int first_bx = bx_ - delayBX;
   int second_bx = 99;
   for (int i = delayBX; i >= 0; i--) {
     if (counter.at(i) >= 2) { // If 2 or more stubs are delayed by i BX or more
       second_bx = bx_ - i; // if i == delayBX, analyze immediately
       break;
     }
   }
   if (not(second_bx != 99))
     { edm::LogError("L1T") << "second_bx = " << second_bx; return; }
 
   // ___________________________________________________________________________
   // Output
 
   track.set_first_bx  ( first_bx );
   track.set_second_bx ( second_bx );
 }

void AngleCalculation::configure	(	int	verbose,
		int	endcap,
		int	sector,
		int	bx,
		int	bxWindow,
		int	thetaWindow,
		int	thetaWindowRPC,
		bool	bugME11Dupes
	)

Definition at line 14 of file AngleCalculation.cc.

References bugME11Dupes_, bx_, bxWindow_, makeMuonMisalignmentScenario::endcap, endcap_, sector_, thetaWindow_, thetaWindowRPC_, verbose, and verbose_.

Referenced by SectorProcessor::process_single_bx().

   {
   verbose_ = verbose;
   endcap_  = endcap;
   sector_  = sector;
   bx_      = bx;
 
   bxWindow_        = bxWindow;
   thetaWindow_     = thetaWindow;
   thetaWindowRPC_  = thetaWindowRPC;
   bugME11Dupes_    = bugME11Dupes;
 }

void AngleCalculation::erase_tracks ( EMTFTrackCollection & tracks ) const

Definition at line 456 of file AngleCalculation.cc.

References emtf::NUM_STATIONS, HiIsolationCommonParameters_cff::track, and x.

Referenced by process().

                                                                      {
   // Erase tracks with rank == 0
   // using erase-remove idiom
   struct {
     typedef EMTFTrack value_type;
     bool operator()(const value_type& x) const {
       return (x.Rank() == 0);
     }
   } rank_zero_pred;
 
   tracks.erase(std::remove_if(tracks.begin(), tracks.end(), rank_zero_pred), tracks.end());
 
   for (const auto& track : tracks) {
     if (not(!track.Hits().empty()))
       { edm::LogError("L1T") << "track.Hits().empty() = " << track.Hits().empty(); return; }
     if (not(track.Hits().size() <= emtf::NUM_STATIONS))
       { edm::LogError("L1T") << "track.Hits().size() = " << track.Hits().size() << ", emtf::NUM_STATIONS= " << emtf::NUM_STATIONS; return; }
   }
 }

void AngleCalculation::process ( emtf::zone_array< EMTFTrackCollection > & zone_tracks ) const

Definition at line 31 of file AngleCalculation.cc.

References calculate_angles(), calculate_bx(), gather_cfg::cout, erase_tracks(), emtf::NUM_ZONES, to_hex(), HiIsolationCommonParameters_cff::track, l1t::tracks, and verbose_.

         {
 
   for (int izone = 0; izone < emtf::NUM_ZONES; ++izone) {
     EMTFTrackCollection& tracks = zone_tracks.at(izone);  // pass by reference
 
     EMTFTrackCollection::iterator tracks_it  = tracks.begin();
     EMTFTrackCollection::iterator tracks_end = tracks.end();
 
     // Calculate deltas
     for (; tracks_it != tracks_end; ++tracks_it) {
       calculate_angles(*tracks_it);
     }
 
     // Erase tracks with rank = 0
     // Erase hits that are not selected as the best phi and theta in each station
     erase_tracks(tracks);
 
     tracks_it  = tracks.begin();
     tracks_end = tracks.end();
 
     // Calculate bx
     // (in the firmware, this happens during best track selection.)
     for (; tracks_it != tracks_end; ++tracks_it) {
       calculate_bx(*tracks_it);
     }
   }  // end loop over zones
 
   if (verbose_ > 0) {  // debug
     for (const auto& tracks : zone_tracks) {
       for (const auto& track : tracks) {
         std::cout << "deltas: z: " << track.Zone()-1 << " pat: " << track.Winner() << " rank: " << to_hex(track.Rank())
             << " delta_ph: " << array_as_string(track.PtLUT().delta_ph)
             << " delta_th: " << array_as_string(track.PtLUT().delta_th)
             << " sign_ph: " << array_as_string(track.PtLUT().sign_ph)
             << " sign_th: " << array_as_string(track.PtLUT().sign_th)
             << " phi: " << track.Phi_fp()
             << " theta: " << track.Theta_fp()
             << " cpat: " << array_as_string(track.PtLUT().cpattern)
             << " v: " << array_as_string(track.PtLUT().bt_vi)
             << " h: " << array_as_string(track.PtLUT().bt_hi)
             << " c: " << array_as_string(track.PtLUT().bt_ci)
             << " s: " << array_as_string(track.PtLUT().bt_si)
             << std::endl;
       }
     }
   }
 
 }