DTTriggerPhase2/src/TrapezoidalGrouping.cc

0001 #include <memory>
0002
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004
0005 #include "L1Trigger/DTTriggerPhase2/interface/TrapezoidalGrouping.h"
0006
0007 using namespace edm;
0008 using namespace std;
0009 using namespace cmsdt;
0010 using namespace dtamgrouping;
0011 // ============================================================================
0012 // Constructors and destructor
0013 // ============================================================================
0014 TrapezoidalGrouping::TrapezoidalGrouping(const ParameterSet &pset, edm::ConsumesCollector &iC)
0015     : MotherGrouping(pset, iC), debug_(pset.getUntrackedParameter<bool>("debug")), currentBaseChannel_(-1) {
0016   // Obtention of parameters
0017   if (debug_)
0018     LogDebug("TrapezoidalGrouping") << "TrapezoidalGrouping: constructor";
0019
0020   // Initialisation of channelIn array
0021   chInDummy_.push_back(DTPrimitive());
0022   for (int lay = 0; lay < NUM_LAYERS; lay++) {
0023     for (int ch = 0; ch < NUM_CH_PER_LAYER; ch++) {
0024       channelIn_[lay][ch] = {chInDummy_};
0025       channelIn_[lay][ch].clear();
0026     }
0027   }
0028 }
0029
0030 TrapezoidalGrouping::~TrapezoidalGrouping() {
0031   if (debug_)
0032     LogDebug("TrapezoidalGrouping") << "TrapezoidalGrouping: destructor";
0033 }
0034
0035 // ============================================================================
0036 // Main methods (initialise, run, finish)
0037 // ============================================================================
0038 void TrapezoidalGrouping::initialise(const edm::EventSetup &iEventSetup) {
0039   if (debug_)
0040     LogDebug("TrapezoidalGrouping") << "TrapezoidalGrouping::initialiase";
0041 }
0042
0043 void TrapezoidalGrouping::run(Event &iEvent,
0044                               const EventSetup &iEventSetup,
0045                               const DTDigiCollection &digis,
0046                               MuonPathPtrs &mpaths) {
0047   //   This function returns the analyzable mpath collection back to the the main function
0048   //   so it can be fitted. This is in fact doing the so-called grouping.
0049   for (int supLayer = 0; supLayer < NUM_SUPERLAYERS; supLayer++) {  // for each SL:
0050     if (debug_)
0051       LogDebug("TrapezoidalGrouping") << "TrapezoidalGrouping::run Reading SL" << supLayer;
0052     setInChannels(&digis, supLayer);
0053
0054     for (auto &hit : all_hits) {
0055       int layer_to_pivot = hit.layerId();
0056       int channel_to_pivot = hit.channelId();
0057       DTPrimitives hits_in_trapezoid;
0058       std::vector<DTPrimitives> hit_mpaths;
0059       std::vector<int> hit_tasks;
0060       for (size_t itask = 0; itask < task_list.size(); itask++) {
0061         // when pivoting over an internal layer, there are two cases
0062         // where the second layer is duplicated
0063         // 12 (0, 5) <-> 14 (0, 7)
0064         // 15 (1, 6) <-> 17 (1, 8)
0065         // we leave it hard-coded here, could be moved somewhere else
0066         if (layer_to_pivot == 1 || layer_to_pivot == 2) {
0067           if (itask == 14 || itask == 17)
0068             continue;
0069         }
0070
0071         auto task = task_list[itask];
0072
0073         std::vector<DTPrimitives> task_mpaths;
0074         std::stack<std::pair<DTPrimitives, int>> mpath_cells_per_task;
0075         mpath_cells_per_task.push(std::make_pair(DTPrimitives({hit}), 0));
0076
0077         while (!mpath_cells_per_task.empty()) {
0078           auto mpath_cells = std::move(mpath_cells_per_task.top());
0079           std::vector<DTPrimitives> tmp_mpaths = {mpath_cells.first};
0080           auto task_index = mpath_cells.second;
0081           auto cell = task[task_index];
0082           auto vertical_shift = trapezoid_vertical_mapping[layer_to_pivot][cell];
0083           auto horizontal_shift = trapezoid_horizontal_mapping[layer_to_pivot][cell];
0084           if (channel_to_pivot + horizontal_shift >= 0 && channel_to_pivot + horizontal_shift < NUM_CH_PER_LAYER) {
0085             tmp_mpaths = group_hits(hit,
0086                                     tmp_mpaths,
0087                                     channelIn_[layer_to_pivot + vertical_shift][channel_to_pivot + horizontal_shift],
0088                                     hits_in_trapezoid);
0089           }
0090           mpath_cells_per_task.pop();
0091           for (const auto &tmp_mpath : tmp_mpaths) {
0092             mpath_cells_per_task.push(std::make_pair(tmp_mpath, task_index + 1));
0093           }
0094           while (!mpath_cells_per_task.empty()) {
0095             if (mpath_cells_per_task.top().second == (int)task.size()) {
0096               task_mpaths.push_back(mpath_cells_per_task.top().first);
0097               mpath_cells_per_task.pop();
0098             } else
0099               break;
0100           }
0101         }
0102         for (auto &task_mpath : task_mpaths) {
0103           hit_mpaths.push_back(task_mpath);
0104           hit_tasks.push_back(itask);
0105         }
0106       }
0107       if (hits_in_trapezoid.size() <= PATHFINDER_INPUT_HITS_LIMIT) {
0108         for (size_t ipath = 0; ipath < hit_mpaths.size(); ipath++) {
0109           auto ptrPrimitive = hit_mpaths[ipath];
0110           auto itask = hit_tasks[ipath];
0111
0112           // In any case, if we have less than 3 hits, we don't output the mpath
0113           if (ptrPrimitive.size() < 3)
0114             continue;
0115
0116           // check if the task has a missing layer associated
0117           // if it does, we add a dummy hit in the missing layer
0118           // if it does not, we check that we actually have 4 present hits;
0119           // if not, we skip the mpath.
0120           if (MISSING_LAYER_LAYOUTS_PER_TASK[layer_to_pivot][itask] != -1) {
0121             auto dtpAux = DTPrimitive();
0122             dtpAux.setTDCTimeStamp(-1);
0123             dtpAux.setChannelId(-1);
0124             dtpAux.setLayerId(MISSING_LAYER_LAYOUTS_PER_TASK[layer_to_pivot][itask]);  //  L=0,1,2,3
0125             dtpAux.setSuperLayerId(hit.superLayerId());
0126             dtpAux.setCameraId(-1);
0127             ptrPrimitive.push_back(dtpAux);
0128           } else {  // we have no missing hits, it must be a 4-hit TP.
0129             if (ptrPrimitive.size() < 4)
0130               continue;
0131           }
0132
0133           // sort the hits by layer, so they are included ordered in the MuonPath object
0134           std::stable_sort(ptrPrimitive.begin(), ptrPrimitive.end(), hitLayerSort);
0135
0136           auto ptrMuonPath = std::make_shared<MuonPath>(ptrPrimitive);
0137           ptrMuonPath->setCellHorizontalLayout(CELL_HORIZONTAL_LAYOUTS_PER_TASK[layer_to_pivot][itask]);
0138           ptrMuonPath->setMissingLayer(MISSING_LAYER_LAYOUTS_PER_TASK[layer_to_pivot][itask]);
0139           mpaths.push_back(std::move(ptrMuonPath));
0140         }
0141       }
0142     }
0143   }
0144   if (debug_)
0145     LogDebug("TrapezoidalGrouping") << "[TrapezoidalGrouping::run] end";
0146 }
0147
0148 void TrapezoidalGrouping::finish() { return; };
0149
0150 // ============================================================================
0151 // Other methods
0152 // ============================================================================
0153
0154 void TrapezoidalGrouping::setInChannels(const DTDigiCollection *digis, int sl) {
0155   //   before setting channels we need to clear
0156   for (int lay = 0; lay < NUM_LAYERS; lay++) {
0157     for (int ch = 0; ch < NUM_CH_PER_LAYER; ch++) {
0158       channelIn_[lay][ch].clear();
0159     }
0160   }
0161   all_hits.clear();
0162
0163   // now fill with those primitives that make sense:
0164   for (const auto &dtLayerId_It : *digis) {
0165     const DTLayerId dtLId = dtLayerId_It.first;
0166
0167     if (dtLId.superlayer() != sl + 1)
0168       continue;  //skip digis not in SL...
0169
0170     for (DTDigiCollection::const_iterator digiIt = (dtLayerId_It.second).first; digiIt != (dtLayerId_It.second).second;
0171          ++digiIt) {
0172       int layer = dtLId.layer() - 1;
0173       int wire = (*digiIt).wire() - 1;
0174       int digiTIME = (*digiIt).time();
0175       int digiTIMEPhase2 = digiTIME;
0176
0177       if (debug_)
0178         LogDebug("TrapezoidalGrouping") << "[TrapezoidalGrouping::setInChannels] SL" << sl << " L" << layer << " : "
0179                                         << wire << " " << digiTIMEPhase2;
0180       auto dtpAux = DTPrimitive();
0181       dtpAux.setTDCTimeStamp(digiTIMEPhase2);
0182       dtpAux.setChannelId(wire);
0183       dtpAux.setLayerId(layer);    //  L=0,1,2,3
0184       dtpAux.setSuperLayerId(sl);  // SL=0,1,2
0185       dtpAux.setCameraId(dtLId.rawId());
0186       channelIn_[layer][wire].push_back(dtpAux);
0187       all_hits.push_back(dtpAux);
0188     }
0189   }
0190
0191   // sort everything by the time of the hits, so it has the same behaviour as the fw
0192   for (int lay = 0; lay < NUM_LAYERS; lay++) {
0193     for (int ch = 0; ch < NUM_CH_PER_LAYER; ch++) {
0194       std::stable_sort(channelIn_[lay][ch].begin(), channelIn_[lay][ch].end(), hitTimeSort);
0195     }
0196   }
0197   std::stable_sort(all_hits.begin(), all_hits.end(), hitTimeSort);
0198 }
0199
0200 std::vector<DTPrimitives> TrapezoidalGrouping::group_hits(DTPrimitive pivot_hit,
0201                                                           std::vector<DTPrimitives> input_paths,
0202                                                           DTPrimitives hits_per_cell,
0203                                                           DTPrimitives &hits_in_trapezoid) {
0204   std::vector<DTPrimitives> output_paths;
0205   for (auto &hit : hits_per_cell) {
0206     int hit_bx = hit.tdcTimeStamp() / LHC_CLK_FREQ;
0207     int pivot_hit_bx = pivot_hit.tdcTimeStamp() / LHC_CLK_FREQ;
0208     if (hitTimeSort(pivot_hit, hit) || (pivot_hit_bx / BX_PER_FRAME) - (hit_bx / BX_PER_FRAME) > MAX_FRAME_DIF)
0209       continue;
0210     // limit the number of hits in the trapezoid to PATHFINDER_INPUT_HITS_LIMIT
0211     if (std::find(hits_in_trapezoid.begin(), hits_in_trapezoid.end(), hit) == hits_in_trapezoid.end())
0212       hits_in_trapezoid.push_back(hit);
0213
0214     if (hits_in_trapezoid.size() > PATHFINDER_INPUT_HITS_LIMIT) {
0215       std::vector<DTPrimitives> empty_paths;
0216       return empty_paths;
0217     }
0218
0219     for (auto &input_path : input_paths) {
0220       auto tmp_path = input_path;
0221       tmp_path.push_back(hit);
0222       output_paths.push_back(tmp_path);
0223     }
0224   }
0225   if (output_paths.empty())
0226     return input_paths;
0227   else
0228     return output_paths;
0229 }