TrackTrigger/src/TTStubAlgorithm_official.cc

0001 /*! \brief   Implementation of methods of TTStubAlgorithm_official
0002  *  \details Here, in the source file, the methods which do depend
0003  *           on the specific type <T> that can fit the template.
0004  *
0005  *  \author Nicola Pozzobon
0006  *  \author Sebastien Viret
0007  *  \date   2013, Jul 18
0008  *
0009  */
0010
0011 #include "L1Trigger/TrackTrigger/interface/TTStubAlgorithm_official.h"
0012
0013 /// Matching operations
0014 template <>
0015 void TTStubAlgorithm_official<Ref_Phase2TrackerDigi_>::PatternHitCorrelation(
0016     bool& aConfirmation,
0017     int& aDisplacement,
0018     int& anOffset,
0019     float& anHardBend,
0020     const TTStub<Ref_Phase2TrackerDigi_>& aTTStub) const {
0021   /// Calculate average coordinates col/row for inner/outer Cluster
0022   // These are already corrected for being at the center of each pixel
0023   MeasurementPoint mp0 = aTTStub.clusterRef(0)->findAverageLocalCoordinates();
0024   MeasurementPoint mp1 = aTTStub.clusterRef(1)->findAverageLocalCoordinates();
0025
0026   bool isPS = aTTStub.moduleTypePS();  // get it from the stub now
0027
0028   /// Get the module position in global coordinates
0029   // TODO temporary: should use a method from the topology
0030   DetId stDetId(aTTStub.getDetId());
0031   const GeomDetUnit* det0 = theTrackerGeom_->idToDetUnit(stDetId + 1);
0032   const GeomDetUnit* det1 = theTrackerGeom_->idToDetUnit(stDetId + 2);
0033
0034   /// Find pixel pitch and topology related information
0035   const PixelGeomDetUnit* pix0 = dynamic_cast<const PixelGeomDetUnit*>(det0);
0036   const PixelGeomDetUnit* pix1 = dynamic_cast<const PixelGeomDetUnit*>(det1);
0037   const PixelTopology* top0 = dynamic_cast<const PixelTopology*>(&(pix0->specificTopology()));
0038   const PixelTopology* top1 = dynamic_cast<const PixelTopology*>(&(pix1->specificTopology()));
0039   std::pair<float, float> pitch0 = top0->pitch();
0040   std::pair<float, float> pitch1 = top1->pitch();
0041
0042   /// Stop if the clusters are not in the same z-segment
0043   int cols0 = top0->ncolumns();
0044   int cols1 = top1->ncolumns();
0045   int ratio = cols0 / cols1;  /// This assumes the ratio is integer!
0046   int segment0 = floor(mp0.y() / ratio);
0047
0048   //  if ( ratio == 1 ) /// 2S Modules
0049   if (!isPS) {
0050     if (mPerformZMatching2S && (segment0 != floor(mp1.y())))
0051       return;
0052   } else  /// PS Modules
0053   {
0054     if (mPerformZMatchingPS && (segment0 != floor(mp1.y())))
0055       return;
0056   }
0057
0058   /// Get the Stack radius and z and displacements
0059   double R0 = det0->position().perp();
0060   double R1 = det1->position().perp();
0061   double Z0 = det0->position().z();
0062   double Z1 = det1->position().z();
0063
0064   double DR = R1 - R0;
0065   double DZ = Z1 - Z0;
0066
0067   double alpha = atan2(DR, DZ);
0068   double delta = sqrt(DR * DR + DZ * DZ) / (R0 * sin(alpha) + Z0 * cos(alpha));
0069
0070   int window = 0;
0071
0072   /// Scale factor is already present in
0073   /// double mPtScalingFactor = (floor(mMagneticFieldStrength*10.0 + 0.5))/10.0*0.0015/mPtThreshold;
0074   /// hence the formula iis something like
0075   /// displacement < Delta * 1 / sqrt( ( 1/(mPtScalingFactor*R) )** 2 - 1 )
0076
0077   /// POSITION IN TERMS OF PITCH MULTIPLES:
0078   ///       0 1 2 3 4 5 5 6 8 9 ...
0079   /// COORD: 0 1 2 3 4 5 6 7 8 9 ...
0080   /// OUT   | | | | | |x| | | | | | | | | |
0081   ///
0082   /// IN    | | | |x|x| | | | | | | | | | |
0083   ///             THIS is 3.5 (COORD) and 4.0 (POS)
0084   /// 1) disp is the difference between average row coordinates
0085   ///    in inner and outer stack member, in terms of outer member pitch
0086   ///    (in case they are the same, this is just a plain coordinate difference)
0087   double dispD = 2 * (mp1.x() - mp0.x()) * (pitch0.first / pitch1.first);  /// In HALF-STRIP units!
0088   int dispI = ((dispD > 0) - (dispD < 0)) * floor(std::abs(dispD));        /// In HALF-STRIP units!
0089   /// 2) offset is the projection with a straight line of the innermost
0090   ///    hit towards the ourermost stack member, still in terms of outer member pitch
0091   ///    NOTE: in terms of coordinates, the center of the module is at NROWS/2-0.5 to
0092   ///    be consistent with the definition given above
0093
0094   /// In HALF-STRIP units!
0095   double offsetD = 2 * delta * (mp0.x() - (top0->nrows() / 2 - 0.5)) * (pitch0.first / pitch1.first);
0096   int offsetI = ((offsetD > 0) - (offsetD < 0)) * floor(std::abs(offsetD));  /// In HALF-STRIP units!
0097
0098   if (stDetId.subdetId() == StripSubdetector::TOB) {
0099     int layer = theTrackerTopo_->layer(stDetId);
0100     int ladder = theTrackerTopo_->tobRod(stDetId);
0101     int type = 2 * theTrackerTopo_->tobSide(stDetId) - 3;  // -1 for tilted-, 1 for tilted+, 3 for flat
0102     double corr = 0;
0103
0104     if (type < 3)  // Only for tilted modules
0105     {
0106       corr = (barrelNTilt.at(layer) + 1) / 2.;
0107       // Corrected ring number, bet 0 and barrelNTilt.at(layer), in ascending |z|
0108       ladder = corr - (corr - ladder) * type;
0109       window = 2 * (tiltedCut.at(layer)).at(ladder);
0110     } else  // Classis barrel window otherwise
0111     {
0112       window = 2 * barrelCut.at(layer);
0113     }
0114
0115   } else if (stDetId.subdetId() == StripSubdetector::TID) {
0116     window = 2 * (ringCut.at(theTrackerTopo_->tidWheel(stDetId))).at(theTrackerTopo_->tidRing(stDetId));
0117   }
0118
0119   /// Accept the stub if the post-offset correction displacement is smaller than the half-window
0120   if (std::abs(dispI - offsetI) <= window)  /// In HALF-STRIP units!
0121   {
0122     aConfirmation = true;
0123     aDisplacement = dispI;                                                     /// In HALF-STRIP units!
0124     anOffset = offsetI;                                                        /// In HALF-STRIP units!
0125     anHardBend = this->degradeBend(isPS, window, (aDisplacement - anOffset));  // In strips units
0126   }  /// End of stub is accepted
0127 }
0128
0129 //--- Does the actual work of degrading the bend. (based on I.Tomalin's code)
0130 template <>
0131 float TTStubAlgorithm_official<Ref_Phase2TrackerDigi_>::degradeBend(bool psModule, int window, int bend) const {
0132   // Number of bits used to encoded bend output by FE electronics.
0133   const unsigned int bitsPS_ = 3;
0134   const unsigned int bits2S_ = 4;
0135
0136   // Number of degraded bend values should correspond to 3 bits (PS modules) or 4 bits (2S modules),
0137   // so measuring everything in half-strip units, max integer "window" size that can be encoded without
0138   // compression given by 2*window+1 <= pow(2,B), where B is number of bits.
0139   // Hence no compression required if window cut is abs(b) <= 3 (PS) or 7 (2S). Must introduce one merge for
0140   // each 1 unit increase in "window" beyond this.
0141
0142   // Bend is measured with granularity of 0.5 strips.
0143   // Convert it to integer measured in half-strip units for this calculation!
0144
0145   float degradedB;
0146   unsigned int numBends = 2 * window + 1;
0147   unsigned int numAllowed = (psModule) ? pow(2, bitsPS_) : pow(2, bits2S_);
0148
0149   // Existance of bend = 0 means can only use an odd number of groups.
0150   numAllowed -= 1;  // NumAllowed can be only based on 3 or 4 bits encoded bends, so 7 or 15 possible values
0151   if (numBends <= numAllowed) {
0152     // Can output uncompressed bend info. (So if window is lower or equal than 1.5 in PS, and 3.5 in 2S
0153     degradedB = static_cast<double>(bend);
0154   } else  // all other cases, need to compress
0155   {
0156     unsigned int inSmallGroup = numBends / numAllowed;
0157     unsigned int numLargeGroups = numBends % numAllowed;
0158     unsigned int inLargeGroup = inSmallGroup + 1;
0159     unsigned int numSmallGroups = numAllowed - numLargeGroups;
0160
0161     std::vector<unsigned int> groups;
0162
0163     // At the end we have
0164     //
0165     // numBends=inSmallGroup*numSmallGroups+inLargeGroup*numLargeGroups
0166     // and
0167     // numAllowed= numSmallGroups+numLargeGroups;
0168     //
0169     // Then you alternate large-small-large-small....large. In the middle, you
0170     // put either large or small, depending if group size is odd or not
0171
0172     for (unsigned int i = 0; i < numLargeGroups / 2; i++)
0173       groups.push_back(inLargeGroup);
0174     for (unsigned int i = 0; i < numSmallGroups / 2; i++)
0175       groups.push_back(inSmallGroup);
0176
0177     // Only one of numLargeGroups & numSmallGroups can be odd, since numAllowed is odd.
0178     // And whichever one is odd is associated to a group with an odd number of elements since numBends is odd,
0179     if (numLargeGroups % 2 == 1 && inLargeGroup % 2 == 1) {
0180       groups.push_back(inLargeGroup);
0181     } else if (numSmallGroups % 2 == 1 && inSmallGroup % 2 == 1) {
0182       groups.push_back(inSmallGroup);
0183     } else {
0184       throw cms::Exception("DegradeBend: logic error with odd numbers");
0185     }
0186
0187     for (unsigned int i = 0; i < numSmallGroups / 2; i++)
0188       groups.push_back(inSmallGroup);
0189     for (unsigned int i = 0; i < numLargeGroups / 2; i++)
0190       groups.push_back(inLargeGroup);
0191
0192     degradedB = 999;
0193     int iUp = -static_cast<int>(window) - 1;  // Start with the minimal possible bend -1
0194     int iDown;
0195
0196     for (unsigned int& inGroup : groups) {
0197       iUp += inGroup;
0198       iDown = iUp - inGroup + 1;
0199       if (bend <= iUp && bend >= iDown) {
0200         degradedB = 0.5 * (iUp + iDown);
0201       }
0202     }
0203     if (degradedB == 999)
0204       throw cms::Exception(
0205           "DegradeStubResolution: error in the group creation, method has been called with wrong inputs");
0206   }
0207
0208   // This is degraded bend in full strip units (neglecting bend sign).
0209   return static_cast<float>(degradedB) / 2.;
0210 }
0211
0212 //--- Check for mistakes.