Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-10-14 01:44:20

0001 //
0002 //  Created by J.Li on 1/23/21.
0003 //
0004 
0005 #include "L1Trigger/TrackTrigger/interface/HitPatternHelper.h"
0006 #include "L1Trigger/TrackFindingTMTT/interface/KFbase.h"
0007 #include "L1Trigger/TrackFindingTMTT/interface/TrackerModule.h"
0008 
0009 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0010 
0011 #include <algorithm>
0012 #include <cmath>
0013 
0014 namespace hph {
0015 
0016   Setup::Setup(const edm::ParameterSet& iConfig, const tt::Setup& setupTT)
0017       : iConfig_(iConfig),
0018         setupTT_(setupTT),
0019         layerIds_(),
0020         layermap_(),
0021         nEtaRegions_(tmtt::KFbase::nEta_ / 2),
0022         nKalmanLayers_(tmtt::KFbase::nKFlayer_) {
0023     for (const tt::SensorModule& sm : setupTT_.sensorModules()) {
0024       layerIds_.push_back(std::make_pair(sm.layerId(), sm.barrel()));
0025     }
0026     sort(layerIds_.begin(), layerIds_.end(), smallerID);
0027     layerIds_.erase(unique(layerIds_.begin(), layerIds_.end(), equalID), layerIds_.end());  //Keep only unique layerIds
0028     // Converting tmtt::KFbase::layerMap_ to a format that is acceptatble by HitPatternHelper
0029     for (int i = 0; i < nEtaRegions_; i++) {
0030       for (int j = 0; j < (int)layerIds_.size(); j++) {
0031         int layer;
0032         if (layerIds_[j].second) {
0033           layer = tmtt::KFbase::layerMap_[i][tmtt::TrackerModule::calcLayerIdReduced(layerIds_[j].first)].first;
0034         } else {
0035           layer = tmtt::KFbase::layerMap_[i][tmtt::TrackerModule::calcLayerIdReduced(layerIds_[j].first)].second;
0036         }
0037         if (layer < nKalmanLayers_) {
0038           layermap_[i][layer].push_back(layerIds_[j].first);
0039         }
0040       }
0041     }
0042   }
0043 
0044   HitPatternHelper::HitPatternHelper(const Setup* setup, int hitpattern, double cot, double z0)
0045       : hitpattern_(hitpattern),
0046         numExpLayer_(0),
0047         numMissingLayer_(0),
0048         numMissingPS_(0),
0049         numMissing2S_(0),
0050         numPS_(0),
0051         num2S_(0),
0052         numMissingInterior1_(0),
0053         numMissingInterior2_(0),
0054         cot_(cot),
0055         z0_(z0),
0056         setup_(setup),
0057         binary_(11, 0),
0058         hphDebug_(setup_->hphDebug()),
0059         useNewKF_(setup_->useNewKF()),
0060         chosenRofZ_(setup_->chosenRofZ()),
0061         deltaTanL_(setup_->deltaTanL()),
0062         etaRegions_(setup_->etaRegions()),
0063         nKalmanLayers_(setup_->nKalmanLayers()),
0064         layermap_(setup_->layermap()) {
0065     //Calculating eta sector based on cot and z0
0066     float kfzRef = z0_ + chosenRofZ_ * cot_;
0067     int kf_eta_reg = 0;
0068     for (int iEtaSec = 1; iEtaSec < ((int)etaRegions_.size() - 1); iEtaSec++) {  // Doesn't apply eta < 2.4 cut.
0069       float etaMax = etaRegions_[iEtaSec];
0070       float zRefMax = chosenRofZ_ / tan(2. * atan(exp(-etaMax)));
0071       if (kfzRef > zRefMax) {
0072         kf_eta_reg = iEtaSec;
0073       }
0074     }
0075     etaSector_ = kf_eta_reg;
0076     if (kf_eta_reg < ((int)etaRegions_.size() - 1) / 2) {
0077       kf_eta_reg = ((int)etaRegions_.size() - 1) / 2 - 1 - kf_eta_reg;
0078     } else {
0079       kf_eta_reg = kf_eta_reg - (int)(etaRegions_.size() - 1) / 2;
0080     }
0081     //Looping over sensor modules to make predictions on which layers particles are expected to hit
0082     for (const tt::SensorModule& sm : setup_->sensorModules()) {
0083       double d = (z0_ - sm.z() + sm.r() * cot_) / (sm.cosTilt() - sm.sinTilt() * cot_);
0084       double d_p =
0085           (z0_ - sm.z() + sm.r() * (cot_ + deltaTanL_ / 2)) / (sm.cosTilt() - sm.sinTilt() * (cot_ + deltaTanL_ / 2));
0086       double d_m =
0087           (z0_ - sm.z() + sm.r() * (cot_ - deltaTanL_ / 2)) / (sm.cosTilt() - sm.sinTilt() * (cot_ - deltaTanL_ / 2));
0088       if (useNewKF_ &&
0089           (abs(d_p) < sm.numColumns() * sm.pitchCol() / 2. || abs(d_m) < sm.numColumns() * sm.pitchCol() / 2.)) {
0090         layers_.push_back(sm);
0091       }
0092       if (!useNewKF_ && abs(d) < sm.numColumns() * sm.pitchCol() / 2.) {
0093         layers_.push_back(sm);
0094       }
0095     }
0096     //layers_ constains all the sensor modules that particles are expected to hit
0097     sort(layers_.begin(), layers_.end(), smallerID);
0098     layers_.erase(unique(layers_.begin(), layers_.end(), equalID), layers_.end());  //Keep only one sensor per layer
0099 
0100     numExpLayer_ = layers_.size();
0101 
0102     int nbits = floor(log2(hitpattern_)) + 1;
0103     int lay_i = 0;
0104     bool seq = false;
0105     for (int i = 0; i < nbits; i++) {
0106       lay_i = ((1 << i) & hitpattern_) >> i;  //0 or 1 in ith bit (right to left)
0107 
0108       if (lay_i && !seq)
0109         seq = true;  //sequence starts when first 1 found
0110       if (!lay_i && seq) {
0111         numMissingInterior1_++;  //This is the same as the "tmp_trk_nlaymiss_interior" calculated in Trackquality.cc
0112       }
0113       if (!lay_i) {
0114         bool realhit = false;
0115         if (layermap_[kf_eta_reg][i].empty())
0116           continue;
0117         for (int j : layermap_[kf_eta_reg][i]) {
0118           int k = findLayer(j);
0119           if (k > 0)
0120             realhit = true;
0121         }
0122         if (realhit)
0123           numMissingInterior2_++;  //This variable doesn't make sense for new KF because it uses the layermap from Old KF
0124       }
0125     }
0126 
0127     if (hphDebug_) {
0128       if (useNewKF_) {
0129         edm::LogVerbatim("TrackTriggerHPH") << "Running with New KF";
0130       } else {
0131         edm::LogVerbatim("TrackTriggerHPH") << "Running with Old KF";
0132       }
0133       edm::LogVerbatim("TrackTriggerHPH") << "======================================================";
0134       edm::LogVerbatim("TrackTriggerHPH")
0135           << "Looking at hitpattern " << std::bitset<7>(hitpattern_) << "; Looping over KF layers:";
0136     }
0137 
0138     if (useNewKF_) {
0139       //New KF uses sensor modules to determine the hitmask already
0140       for (int i = 0; i < numExpLayer_; i++) {
0141         if (hphDebug_) {
0142           edm::LogVerbatim("TrackTriggerHPH") << "--------------------------";
0143           edm::LogVerbatim("TrackTriggerHPH") << "Looking at KF layer " << i;
0144           if (layers_[i].layerId() < 10) {
0145             edm::LogVerbatim("TrackTriggerHPH") << "KF expects L" << layers_[i].layerId();
0146           } else {
0147             edm::LogVerbatim("TrackTriggerHPH") << "KF expects D" << layers_[i].layerId() - 10;
0148           }
0149         }
0150 
0151         if (((1 << i) & hitpattern_) >> i) {
0152           if (hphDebug_) {
0153             edm::LogVerbatim("TrackTriggerHPH") << "Layer found in hitpattern";
0154           }
0155 
0156           binary_[ReducedId(layers_[i].layerId())] = 1;
0157           if (layers_[i].psModule()) {
0158             numPS_++;
0159           } else {
0160             num2S_++;
0161           }
0162         } else {
0163           if (hphDebug_) {
0164             edm::LogVerbatim("TrackTriggerHPH") << "Layer missing in hitpattern";
0165           }
0166 
0167           if (layers_[i].psModule()) {
0168             numMissingPS_++;
0169           } else {
0170             numMissing2S_++;
0171           }
0172         }
0173       }
0174 
0175     } else {
0176       //Old KF uses the hard coded layermap to determien hitmask
0177       for (int i = 0; i < nKalmanLayers_; i++) {  //Loop over each digit of hitpattern
0178 
0179         if (hphDebug_) {
0180           edm::LogVerbatim("TrackTriggerHPH") << "--------------------------";
0181           edm::LogVerbatim("TrackTriggerHPH") << "Looking at KF layer " << i;
0182         }
0183 
0184         if (layermap_[kf_eta_reg][i].empty()) {
0185           if (hphDebug_) {
0186             edm::LogVerbatim("TrackTriggerHPH") << "KF does not expect this layer";
0187           }
0188 
0189           continue;
0190         }
0191 
0192         for (int j :
0193              layermap_[kf_eta_reg][i]) {  //Find out which layer the Old KF is dealing with when hitpattern is encoded
0194 
0195           if (hphDebug_) {
0196             if (j < 10) {
0197               edm::LogVerbatim("TrackTriggerHPH") << "KF expects L" << j;
0198             } else {
0199               edm::LogVerbatim("TrackTriggerHPH") << "KF expects D" << j - 10;
0200             }
0201           }
0202 
0203           int k = findLayer(j);
0204           if (k < 0) {
0205             //k<0 means even though layer j is predicted by Old KF, this prediction is rejected because it contradicts
0206             if (hphDebug_) {  //a more accurate prediction made with the help of information from sensor modules
0207               edm::LogVerbatim("TrackTriggerHPH") << "Rejected by sensor modules";
0208             }
0209 
0210             continue;
0211           }
0212 
0213           if (hphDebug_) {
0214             edm::LogVerbatim("TrackTriggerHPH") << "Confirmed by sensor modules";
0215           }
0216           //prediction is accepted
0217           if (((1 << i) & hitpattern_) >> i) {
0218             if (hphDebug_) {
0219               edm::LogVerbatim("TrackTriggerHPH") << "Layer found in hitpattern";
0220             }
0221 
0222             binary_[ReducedId(j)] = 1;
0223             if (layers_[k].psModule()) {
0224               numPS_++;
0225             } else {
0226               num2S_++;
0227             }
0228           } else {
0229             if (hphDebug_) {
0230               edm::LogVerbatim("TrackTriggerHPH") << "Layer missing in hitpattern";
0231             }
0232 
0233             if (layers_[k].psModule()) {
0234               numMissingPS_++;
0235             } else {
0236               numMissing2S_++;
0237             }
0238           }
0239         }
0240       }
0241     }
0242 
0243     if (hphDebug_) {
0244       edm::LogVerbatim("TrackTriggerHPH") << "------------------------------";
0245       edm::LogVerbatim("TrackTriggerHPH") << "numPS = " << numPS_ << ", num2S = " << num2S_
0246                                           << ", missingPS = " << numMissingPS_ << ", missing2S = " << numMissing2S_;
0247       edm::LogVerbatim("TrackTriggerHPH") << "======================================================";
0248     }
0249   }
0250 
0251   int HitPatternHelper::ReducedId(int layerId) {
0252     if (hphDebug_ && (layerId > 15 || layerId < 1)) {
0253       edm::LogVerbatim("TrackTriggerHPH") << "Warning: invalid layer id !";
0254     }
0255     if (layerId <= 6) {
0256       layerId = layerId - 1;
0257       return layerId;
0258     } else {
0259       layerId = layerId - 5;
0260       return layerId;
0261     }
0262   };
0263 
0264   int HitPatternHelper::findLayer(int layerId) {
0265     for (int i = 0; i < (int)layers_.size(); i++) {
0266       if (layerId == (int)layers_[i].layerId()) {
0267         return i;
0268       }
0269     }
0270     return -1;
0271   }
0272 
0273 }  // namespace hph