Local/src/FastLineRecognition.cc

0001 /****************************************************************************
0002 *
0003 * This is a part of TOTEM offline software.
0004 * Authors:
0005 *   Jan Kašpar (jan.kaspar@gmail.com)
0006 *
0007 ****************************************************************************/
0008
0009 #include "RecoPPS/Local/interface/FastLineRecognition.h"
0010
0011 #include "DataFormats/CTPPSReco/interface/TotemRPRecHit.h"
0012
0013 #include "Geometry/VeryForwardGeometryBuilder/interface/CTPPSGeometry.h"
0014
0015 #include <map>
0016 #include <cmath>
0017 #include <cstdio>
0018 #include <algorithm>
0019
0020 //#define CTPPS_DEBUG 1
0021
0022 using namespace std;
0023 using namespace edm;
0024
0025 //----------------------------------------------------------------------------------------------------
0026
0027 const double FastLineRecognition::sigma0 = 66E-3 / sqrt(12.);
0028
0029 //----------------------------------------------------------------------------------------------------
0030
0031 void FastLineRecognition::Cluster::add(const Point *p1, const Point *p2, double a, double b, double w) {
0032   // which points to be added to contents?
0033   bool add1 = true, add2 = true;
0034   for (vector<const Point *>::const_iterator it = contents.begin(); it != contents.end() && (add1 || add2); ++it) {
0035     if ((*it)->hit == p1->hit)
0036       add1 = false;
0037
0038     if ((*it)->hit == p2->hit)
0039       add2 = false;
0040   }
0041
0042   // add the points
0043   if (add1)
0044     contents.push_back(p1);
0045   if (add2)
0046     contents.push_back(p2);
0047
0048   // update sums
0049   Saw += a * w;
0050   Sbw += b * w;
0051   Sw += w;
0052   S1 += 1.;
0053 }
0054
0055 //----------------------------------------------------------------------------------------------------
0056 //----------------------------------------------------------------------------------------------------
0057
0058 FastLineRecognition::FastLineRecognition(double cw_a, double cw_b)
0059     : chw_a(cw_a / 2.), chw_b(cw_b / 2.), geometry(nullptr) {}
0060
0061 //----------------------------------------------------------------------------------------------------
0062
0063 FastLineRecognition::~FastLineRecognition() {}
0064
0065 //----------------------------------------------------------------------------------------------------
0066
0067 FastLineRecognition::GeomData FastLineRecognition::getGeomData(unsigned int id) {
0068   // result already buffered?
0069   map<unsigned int, GeomData>::iterator it = geometryMap.find(id);
0070   if (it != geometryMap.end())
0071     return it->second;
0072
0073   // calculate it
0074   const auto &d = geometry->localToGlobalDirection(id, CTPPSGeometry::Vector(0., 1., 0.));
0075   DetGeomDesc::Translation c = geometry->sensor(TotemRPDetId(id))->translation();
0076   GeomData gd;
0077   gd.z = c.z();
0078   gd.s = d.x() * c.x() + d.y() * c.y();
0079
0080   geometryMap[id] = gd;
0081
0082   return gd;
0083 }
0084
0085 //----------------------------------------------------------------------------------------------------
0086
0087 void FastLineRecognition::getPatterns(const DetSetVector<TotemRPRecHit> &input,
0088                                       double z0,
0089                                       double threshold,
0090                                       DetSet<TotemRPUVPattern> &patterns) {
0091   // build collection of points in the global coordinate system
0092   std::vector<Point> points;
0093   for (auto &ds : input) {
0094     unsigned int detId = ds.detId();
0095
0096     for (auto &h : ds) {
0097       const TotemRPRecHit *hit = &h;
0098       const GeomData &gd = getGeomData(detId);
0099
0100       double p = hit->position() + gd.s;
0101       double z = gd.z - z0;
0102       double w = sigma0 / hit->sigma();
0103
0104       points.push_back(Point(detId, hit, p, z, w));
0105     }
0106   }
0107
0108 #if CTPPS_DEBUG > 0
0109   printf(">> FastLineRecognition::getPatterns(z0 = %E)\n", z0);
0110   printf(">>>>>>>>>>>>>>>>>>>\n");
0111 #endif
0112
0113   // reset output
0114   patterns.clear();
0115
0116   Cluster c;
0117   while (getOneLine(points, threshold, c)) {
0118     // convert cluster to pattern and save it
0119     TotemRPUVPattern pattern;
0120     pattern.setA(c.Saw / c.Sw);
0121     pattern.setB(c.Sbw / c.Sw);
0122     pattern.setW(c.weight);
0123
0124 #if CTPPS_DEBUG > 0
0125     printf("\tpoints of the selected cluster: %lu\n", c.contents.size());
0126 #endif
0127
0128     for (auto &pit : c.contents) {
0129 #if CTPPS_DEBUG > 0
0130       printf("\t\t%.1f\n", pit->z);
0131 #endif
0132       pattern.addHit(pit->detId, *(pit->hit));
0133     }
0134
0135     patterns.push_back(pattern);
0136
0137 #if CTPPS_DEBUG > 0
0138     unsigned int u_points_b = 0;
0139     for (vector<Point>::iterator dit = points.begin(); dit != points.end(); ++dit)
0140       if (dit->usable)
0141         u_points_b++;
0142     printf("\tusable points before: %u\n", u_points_b);
0143 #endif
0144
0145     // remove points belonging to the recognized line
0146     for (vector<const Point *>::iterator hit = c.contents.begin(); hit != c.contents.end(); ++hit) {
0147       for (vector<Point>::iterator dit = points.begin(); dit != points.end(); ++dit) {
0148         //printf("\t\t1: %.2f, %p vs. 2: %.2f, %p\n", (*hit)->z, (*hit)->hit, dit->z, dit->hit);
0149         if ((*hit)->hit == dit->hit) {
0150           dit->usable = false;
0151           //points.erase(dit);
0152           break;
0153         }
0154       }
0155     }
0156
0157 #if CTPPS_DEBUG > 0
0158     unsigned int u_points_a = 0;
0159     for (vector<Point>::iterator dit = points.begin(); dit != points.end(); ++dit)
0160       if (dit->usable)
0161         u_points_a++;
0162     printf("\tusable points after: %u\n", u_points_a);
0163 #endif
0164   }
0165
0166 #if CTPPS_DEBUG > 0
0167   printf("patterns at end: %lu\n", patterns.size());
0168   printf("<<<<<<<<<<<<<<<<<<<\n");
0169 #endif
0170 }
0171
0172 //----------------------------------------------------------------------------------------------------
0173
0174 bool FastLineRecognition::getOneLine(const vector<FastLineRecognition::Point> &points,
0175                                      double threshold,
0176                                      FastLineRecognition::Cluster &result) {
0177 #if CTPPS_DEBUG > 0
0178   printf("\tFastLineRecognition::getOneLine\n");
0179 #endif
0180
0181   if (points.size() < 2)
0182     return false;
0183
0184   vector<Cluster> clusters;
0185
0186   // go through all the combinations of measured points
0187   for (vector<Point>::const_iterator it1 = points.begin(); it1 != points.end(); ++it1) {
0188     if (!it1->usable)
0189       continue;
0190
0191     for (vector<Point>::const_iterator it2 = it1; it2 != points.end(); ++it2) {
0192       if (!it2->usable)
0193         continue;
0194
0195       const double &z1 = it1->z;
0196       const double &z2 = it2->z;
0197
0198       if (z1 == z2)
0199         continue;
0200
0201       const double &p1 = it1->h;
0202       const double &p2 = it2->h;
0203
0204       const double &w1 = it1->w;
0205       const double &w2 = it2->w;
0206
0207       // calculate intersection
0208       double a = (p2 - p1) / (z2 - z1);
0209       double b = p1 - z1 * a;
0210       double w = w1 + w2;
0211
0212 #if CTPPS_DEBUG > 0
0213       printf("\t\t\tz: 1=%+5.1f, 2=%+5.1f | U/V: 1=%+6.3f, 2=%+6.3f | a=%+6.3f rad, b=%+6.3f mm, w=%.1f\n",
0214              z1,
0215              z2,
0216              p1,
0217              p2,
0218              a,
0219              b,
0220              w);
0221 #endif
0222
0223       // add it to the appropriate cluster
0224       bool newCluster = true;
0225       for (unsigned int k = 0; k < clusters.size(); k++) {
0226         Cluster &c = clusters[k];
0227         if (c.S1 < 1. || c.Sw <= 0.)
0228           continue;
0229
0230 #if CTPPS_DEBUG > 0
0231         if (k < 10)
0232           printf("\t\t\t\ttest cluster %u at a=%+6.3f, b=%+6.3f : %+6.3f, %+6.3f : %i, %i\n",
0233                  k,
0234                  c.Saw / c.Sw,
0235                  c.Sbw / c.Sw,
0236                  chw_a,
0237                  chw_b,
0238                  (std::abs(a - c.Saw / c.Sw) < chw_a),
0239                  (std::abs(b - c.Sbw / c.Sw) < chw_b));
0240 #endif
0241
0242         if ((std::abs(a - c.Saw / c.Sw) < chw_a) && (std::abs(b - c.Sbw / c.Sw) < chw_b)) {
0243           newCluster = false;
0244           clusters[k].add(&(*it1), &(*it2), a, b, w);
0245 #if CTPPS_DEBUG > 0
0246           printf("\t\t\t\t--> cluster %u\n", k);
0247 #endif
0248           break;
0249         }
0250       }
0251
0252       // make new cluster
0253       if (newCluster) {
0254 #if CTPPS_DEBUG > 0
0255         printf("\t\t\t\t--> new cluster %lu\n", clusters.size());
0256 #endif
0257         clusters.push_back(Cluster());
0258         clusters.back().add(&(*it1), &(*it2), a, b, w);
0259       }
0260     }
0261   }
0262
0263 #if CTPPS_DEBUG > 0
0264   printf("\t\tclusters: %lu\n", clusters.size());
0265 #endif
0266
0267   // find the cluster with highest weight
0268   unsigned int mk = 0;
0269   double mw = -1.;
0270   for (unsigned int k = 0; k < clusters.size(); k++) {
0271     double w = 0;
0272     for (vector<const Point *>::iterator it = clusters[k].contents.begin(); it != clusters[k].contents.end(); ++it)
0273       w += (*it)->w;
0274     clusters[k].weight = w;
0275
0276     if (w > mw) {
0277       mw = w;
0278       mk = k;
0279     }
0280   }
0281
0282 #if CTPPS_DEBUG > 0
0283   printf("\t\tmw = %.1f, mk = %u\n", mw, mk);
0284 #endif
0285
0286   // rerturn result
0287   if (mw >= threshold) {
0288     result = clusters[mk];
0289
0290     return true;
0291   } else
0292     return false;
0293 }