TrackFindingTMTT/src/TrkRZfilter.cc

0001 #include "FWCore/Utilities/interface/Exception.h"
0002
0003 #include "L1Trigger/TrackFindingTMTT/interface/TrkRZfilter.h"
0004 #include "L1Trigger/TrackFindingTMTT/interface/Settings.h"
0005 #include "L1Trigger/TrackFindingTMTT/interface/Stub.h"
0006 #include "L1Trigger/TrackFindingTMTT/interface/PrintL1trk.h"
0007
0008 #include <array>
0009
0010 using namespace std;
0011
0012 namespace tmtt {
0013
0014   //=== Initialize configuration parameters, and note eta range covered by sector and phi coordinate of its centre.
0015
0016   TrkRZfilter::TrkRZfilter(const Settings* settings,
0017                            unsigned int iPhiSec,
0018                            unsigned int iEtaReg,
0019                            float etaMinSector,
0020                            float etaMaxSector,
0021                            float phiCentreSector)
0022       :  // Configuration parameters.
0023         settings_(settings),
0024         // Sector number.
0025         iPhiSec_(iPhiSec),
0026         iEtaReg_(iEtaReg) {
0027     // Eta range of sector & phi coord of its centre.
0028     etaMinSector_ = etaMinSector;
0029     etaMaxSector_ = etaMaxSector;
0030     phiCentreSector_ = phiCentreSector;
0031
0032     // Note that no estimate of the r-z helix params yet exists.
0033     rzHelix_set_ = false;  // No valid estimate yet.
0034
0035     // Calculate z coordinate a track would have at given radius if on eta sector boundary.
0036     chosenRofZ_ = settings->chosenRofZ();
0037     zTrkMinSector_ = chosenRofZ_ / tan(2. * atan(exp(-etaMinSector_)));
0038     zTrkMaxSector_ = chosenRofZ_ / tan(2. * atan(exp(-etaMaxSector_)));
0039     unsigned int zbits = settings->zBits();
0040     unsigned int rtbits = settings_->rtBits();
0041     float zrange = settings_->zRange();
0042     float rtRange = settings_->rtRange();
0043     float zMultiplier = pow(2., zbits) / zrange;
0044     float rMultiplier = pow(2., rtbits) / rtRange;
0045
0046     if (settings_->enableDigitize()) {
0047       zTrkMinSector_ = floor(zTrkMinSector_ * zMultiplier) / zMultiplier;
0048       zTrkMaxSector_ = floor(zTrkMaxSector_ * zMultiplier) / zMultiplier;
0049       chosenRofZ_ = floor(chosenRofZ_ * rMultiplier) / rMultiplier;
0050     }
0051
0052     // Assumed length of beam-spot in z.
0053     beamWindowZ_ = settings->beamWindowZ();
0054
0055     // Name of r-z track filter algorithm to run.
0056     rzFilterName_ = settings->rzFilterName();
0057
0058     // --- Options for Seed filter.
0059     // Keep stubs compatible with all possible good seed.
0060     keepAllSeed_ = settings->keepAllSeed();
0061     //Added resolution for a tracklet-like filter algorithm, beyond that estimated from hit resolution.
0062     seedResCut_ = settings->seedResCut();
0063     // Maximum number of seed combinations to bother checking per track candidate.
0064     maxSeedCombinations_ = settings->maxSeedCombinations();
0065     // Maximum number of seed combinations consistent with sector (z0,eta) constraints to bother checking per track candidate.
0066     maxGoodSeedCombinations_ = settings->maxGoodSeedCombinations();
0067     // Maximum number of seeds that a single stub can be included in.
0068     maxSeedsPerStub_ = settings->maxSeedsPerStub();
0069     // Reject tracks whose estimated rapidity from seed filter is inconsistent range of with eta sector. (Kills some duplicate tracks).
0070     zTrkSectorCheck_ = settings->zTrkSectorCheck();
0071
0072     // For debugging
0073     minNumMatchLayers_ = settings->minNumMatchLayers();
0074   }
0075
0076   //=== Filters track candidates (found by the r-phi Hough transform), removing inconsistent stubs from the tracks,
0077   //=== also killing some of the tracks altogether if they are left with too few stubs.
0078   //=== Also adds an estimate of r-z helix parameters to the selected track objects, if the filters used provide this.
0079
0080   list<L1track3D> TrkRZfilter::filterTracks(const list<L1track2D>& tracks) {
0081     list<L1track3D> filteredTracks;
0082
0083     for (const L1track2D& trkIN : tracks) {
0084       const vector<Stub*>& stubs = trkIN.stubs();  // stubs assigned to track
0085
0086       // Declare this to be worth keeping (for now).
0087       bool trackAccepted = true;
0088
0089       // Note that no estimate of the r-z helix params of this track yet exists.
0090       rzHelix_set_ = false;  // No valid estimate yet.
0091
0092       // Digitize stubs for r-z filter if required.
0093       if (settings_->enableDigitize()) {
0094         for (Stub* s : stubs) {
0095           s->digitize(iPhiSec_, Stub::DigiStage::SF);
0096         }
0097       }
0098
0099       //--- Filter stubs assigned to track, checking they are consistent with requested criteria.
0100
0101       // Get debug printout for specific regions.
0102       bool print = false;
0103
0104       vector<Stub*> filteredStubs = stubs;
0105       if (rzFilterName_ == "SeedFilter") {
0106         filteredStubs = this->seedFilter(filteredStubs, trkIN.qOverPt(), print);
0107       } else {
0108         throw cms::Exception("BadConfig") << "TrkRzFilter: ERROR unknown r-z track filter requested: " << rzFilterName_;
0109       }
0110
0111       // Check if track still has stubs in enough layers after filter.
0112
0113       unsigned int numLayersAfterFilters = Utility::countLayers(settings_, filteredStubs);
0114       if (numLayersAfterFilters < settings_->minFilterLayers())
0115         trackAccepted = false;
0116       //if (numLayersAfterFilters < Utility::numLayerCut(Utility::AlgoStep::SEED, settings_, iPhiSec_, iEtaReg_, std::abs(trkIN.qOverPt()))) trackAccepted = false;
0117
0118       if (trackAccepted) {
0119         // Estimate r-z helix parameters from centre of eta-sector if no estimate provided by r-z filter.
0120         if (!rzHelix_set_)
0121           this->estRZhelix();
0122
0123         pair<float, float> helixRZ(rzHelix_z0_, rzHelix_tanL_);
0124
0125         // Create copy of original track, except now using its filtered stubs, to be added to filteredTrack collection.
0126         filteredTracks.emplace_back(settings_,
0127                                     filteredStubs,
0128                                     trkIN.cellLocationHT(),
0129                                     trkIN.helix2D(),
0130                                     helixRZ,
0131                                     trkIN.iPhiSec(),
0132                                     trkIN.iEtaReg(),
0133                                     trkIN.optoLinkID(),
0134                                     trkIN.mergedHTcell());
0135       }
0136     }
0137
0138     return filteredTracks;
0139   }
0140
0141   //=== Use Seed Filter to produce a filtered collection of stubs on this track candidate that are consistent with a straight line
0142   //=== in r-z using tracklet algo.
0143
0144   vector<Stub*> TrkRZfilter::seedFilter(const std::vector<Stub*>& stubs, float trkQoverPt, bool print) {
0145     unsigned int numLayers;                //Num of Layers in the cell after that filter has been applied
0146     std::vector<Stub*> filtStubs = stubs;  // Copy stubs vector in filtStubs
0147     bool FirstSeed = true;
0148     //Allowed layers for the first & second seeding stubs
0149     constexpr std::array<unsigned int, 8> FirstSeedLayers = {{1, 2, 11, 21, 3, 12, 22, 4}};
0150     constexpr std::array<unsigned int, 10> SecondSeedLayers = {{1, 2, 11, 3, 21, 22, 12, 23, 13, 4}};
0151     set<Stub*> uniqueFilteredStubs;
0152
0153     unsigned int numSeedCombinations = 0;      // Counter for number of seed combinations considered.
0154     unsigned int numGoodSeedCombinations = 0;  // Counter for seed combinations with z0 within beam spot length.
0155     vector<Stub*> filteredStubs;               // Filter Stubs vector to be returned
0156
0157     unsigned int oldNumLay = 0;  //Number of Layers counter, used to keep the seed with more layers
0158
0159     auto orderByLayer = [](const Stub* a, const Stub* b) { return bool(a->layerId() < b->layerId()); };
0160     std::sort(filtStubs.begin(), filtStubs.end(), orderByLayer);
0161
0162     // Loop over stubs in the HT Cell
0163     for (Stub* s0 : filtStubs) {
0164       // Select the first available seeding stub (r<70)
0165       if (s0->psModule() && std::find(std::begin(FirstSeedLayers), std::end(FirstSeedLayers), s0->layerId()) !=
0166                                 std::end(FirstSeedLayers)) {
0167         unsigned int numSeedsPerStub = 0;
0168
0169         for (Stub* s1 : filtStubs) {
0170           if (numGoodSeedCombinations < maxGoodSeedCombinations_ && numSeedCombinations < maxSeedCombinations_ &&
0171               numSeedsPerStub < maxSeedsPerStub_) {
0172             // Select the second seeding stub (r<90)
0173             if (s1->psModule() && s1->layerId() > s0->layerId() &&
0174                 std::find(std::begin(SecondSeedLayers), std::end(SecondSeedLayers), s1->layerId()) !=
0175                     std::end(SecondSeedLayers)) {
0176               numSeedsPerStub++;
0177               numSeedCombinations++;  //Increase filter cycles counter
0178               if (print)
0179                 PrintL1trk() << "s0: "
0180                              << "z: " << s0->z() << ", r: " << s0->r() << ", id:" << s0->layerId() << " ****** s1: "
0181                              << "z: " << s1->z() << ", r: " << s1->r() << ", id:" << s1->layerId();
0182               //double sumSeedDist = 0.;
0183               //double oldSumSeedDist = 1000000.;  //Define variable used to estimate the quality of seeds
0184               vector<Stub*> tempStubs;  //Create a temporary container for stubs
0185               tempStubs.push_back(s0);  //Store the first seeding stub in the temporary container
0186               tempStubs.push_back(s1);  //Store the second seeding stub in the temporary container
0187
0188               // Estimate a value of z at the beam spot using the two seeding stubs
0189               double z0 = s1->z() + (-s1->z() + s0->z()) * s1->r() / (s1->r() - s0->r());
0190               // COMMENTED OUT CODE GIVES OPTION OF ALLOWING FOR UNCERTAINTY. NOT WORTH IT?
0191               //double z0err = s1->sigmaZ() + ( s1->sigmaZ() + s0->sigmaZ() )*s1->r()/std::abs(s1->r()-s0->r()) + std::abs(-s1->z()+s0->z())*(s1->sigmaR()*std::abs(s1->r()-s0->r()) + s1->r()*(s1->sigmaR() + s0->sigmaR()) )/((s1->r()-s0->r())*(s1->r()-s0->r()));
0192               // Estimate a value of z at a chosen Radius using the two seeding stubs
0193               float zTrk = s1->z() + (-s1->z() + s0->z()) * (s1->r() - chosenRofZ_) / (s1->r() - s0->r());
0194               // float zTrkErr = s1->sigmaZ() + ( s1->sigmaZ() + s0->sigmaZ() )*std::abs(s1->r()-chosenRofZ_)/std::abs(s1->r()-s0->r()) + std::abs(-s1->z()+s0->z())*(s1->sigmaR()*std::abs(s1->r()-s0->r()) + std::abs(s1->r()-chosenRofZ_)*(s1->sigmaR() + s0->sigmaR()) )/((s1->r()-s0->r())*(s1->r()-s0->r()));
0195               float leftZtrk = zTrk * std::abs(s1->r() - s0->r());
0196               float rightZmin = zTrkMinSector_ * std::abs(s1->r() - s0->r());
0197               float rightZmax = zTrkMaxSector_ * std::abs(s1->r() - s0->r());
0198
0199               // If z0 is within the beamspot range loop over the other stubs in the cell
0200               //if (std::abs(z0)<=beamWindowZ_+z0err) {
0201               if (std::abs(z0) <= beamWindowZ_) {
0202                 // Check track r-z helix parameters are consistent with it being assigned to current rapidity sector (kills duplicates due to overlapping sectors).
0203                 // if ( (! zTrkSectorCheck_) || (zTrk > zTrkMinSector_ - zTrkErr && zTrk < zTrkMaxSector_ + zTrkErr) ) {
0204                 if ((!zTrkSectorCheck_) || (leftZtrk > rightZmin && leftZtrk < rightZmax)) {
0205                   numGoodSeedCombinations++;
0206                   unsigned int LiD = 0;  //Store the layerId of the stub (KEEP JUST ONE STUB PER LAYER)
0207                   //                                double oldseed = 1000.; //Store the seed value of the current stub (KEEP JUST ONE STUB PER LAYER)
0208
0209                   // Loop over stubs in vector different from the seeding stubs
0210                   for (Stub* s : filtStubs) {
0211                     // if(s!= s0 && s!= s1){
0212                     // Calculate the seed and its tolerance
0213                     double seedDist =
0214                         (s->z() - s1->z()) * (s1->r() - s0->r()) - (s->r() - s1->r()) * (s1->z() - s0->z());
0215                     double seedDistRes = (s->sigmaZ() + s1->sigmaZ()) * std::abs(s1->r() - s0->r()) +
0216                                          (s->sigmaR() + s1->sigmaR()) * std::abs(s1->z() - s0->z()) +
0217                                          (s0->sigmaZ() + s1->sigmaZ()) * std::abs(s->r() - s1->r()) +
0218                                          (s0->sigmaR() + s1->sigmaR()) * std::abs(s->z() - s1->z());
0219                     seedDistRes *= seedResCut_;  // Scale cut by this amount.
0220                     //If seed is lower than the tolerance push back the stub (KEEP JUST ONE STUB PER LAYER, NOT ENABLED BY DEFAULT)
0221                     if (std::abs(seedDist) <= seedDistRes) {
0222                       if (s->layerId() != LiD and s->layerId() != s0->layerId() and s->layerId() != s1->layerId()) {
0223                         tempStubs.push_back(s);
0224                         LiD = s->layerId();
0225                       }
0226                     }
0227                   }
0228                 }
0229               }
0230
0231               numLayers = Utility::countLayers(
0232                   settings_, tempStubs);  // Count the number of layers in the temporary stubs container
0233
0234               // Check if the current seed has more layers then the previous one (Keep the best seed)
0235               if (keepAllSeed_ == false) {
0236                 if (numLayers > oldNumLay) {
0237                   // Check if the current seed has better quality than the previous one
0238                   //if(sumSeedDist < oldSumSeedDist){
0239                   filteredStubs =
0240                       tempStubs;  //Copy the temporary stubs vector in the filteredStubs vector, which will be returned
0241                                   //  oldSumSeedDist = sumSeedDist; //Update value of oldSumSeedDist
0242                   oldNumLay = numLayers;                                      //Update value of oldNumLay
0243                   rzHelix_z0_ = z0;                                           //Store estimated z0
0244                   rzHelix_tanL_ = (s1->z() - s0->z()) / (s1->r() - s0->r());  // Store estimated tanLambda
0245                   rzHelix_set_ = true;
0246                 }
0247                 //}
0248               } else {
0249                 // Check if the current seed satisfies the minimum layers requirement (Keep all seed algorithm)
0250                 if (numLayers >= Utility::numLayerCut(
0251                                      Utility::AlgoStep::SEED, settings_, iPhiSec_, iEtaReg_, std::abs(trkQoverPt))) {
0252                   uniqueFilteredStubs.insert(tempStubs.begin(), tempStubs.end());  //Insert the uniqueStub set
0253
0254                   // If these are the first seeding stubs store the values of z0 and tanLambda
0255                   if (FirstSeed) {
0256                     FirstSeed = false;
0257                     rzHelix_z0_ = z0;                                           //Store estimated z0
0258                     rzHelix_tanL_ = (s1->z() - s0->z()) / (s1->r() - s0->r());  // Store estimated tanLambda
0259                     rzHelix_set_ = true;
0260                   }
0261                 }
0262               }
0263             }
0264           }
0265         }
0266       }
0267     }
0268
0269     // Copy stubs from the uniqueFilteredStubs set to the filteredStubs vector (Keep all seed algorithm)
0270     if (keepAllSeed_ == true) {
0271       for (Stub* stub : uniqueFilteredStubs) {
0272         filteredStubs.push_back(stub);
0273       }
0274     }
0275
0276     // Note number of seed combinations used for this track.
0277     numSeedCombsPerTrk_.push_back(numSeedCombinations);
0278     numGoodSeedCombsPerTrk_.push_back(numGoodSeedCombinations);
0279
0280     return filteredStubs;  // Return the filteredStubs vector
0281   }
0282
0283   // Estimate r-z helix parameters from centre of eta-sector if no better estimate provided by r-z filter.
0284
0285   void TrkRZfilter::estRZhelix() {
0286     rzHelix_z0_ = 0.;
0287     rzHelix_tanL_ = 0.5 * (1 / tan(2 * atan(exp(-etaMinSector_))) + 1 / tan(2 * atan(exp(-etaMaxSector_))));
0288     rzHelix_set_ = true;
0289   }
0290
0291 }  // namespace tmtt