CSCRecHit/src/CSCSegment.cc

0001 /** \file CSCSegment.cc
0002  *
0003  *  \author Matteo Sani
0004  */
0005
0006 #include <DataFormats/CSCRecHit/interface/CSCSegment.h>
0007 #include <iostream>
0008
0009 namespace {
0010   // Get CSCDetId from one of the rechits, but then remove the layer part so it's a _chamber_ id
0011   inline DetId buildDetId(CSCDetId id) { return CSCDetId(id.endcap(), id.station(), id.ring(), id.chamber(), 0); }
0012
0013 }  // namespace
0014
0015 CSCSegment::CSCSegment(const std::vector<const CSCRecHit2D*>& proto_segment,
0016                        LocalPoint origin,
0017                        LocalVector direction,
0018                        const AlgebraicSymMatrix& errors,
0019                        double chi2)
0020     : RecSegment(buildDetId(proto_segment.front()->cscDetId())),
0021       theOrigin(origin),
0022       theLocalDirection(direction),
0023       theCovMatrix(errors),
0024       theChi2(chi2),
0025       aME11a_duplicate(false) {
0026   for (unsigned int i = 0; i < proto_segment.size(); ++i)
0027     theCSCRecHits.push_back(*proto_segment[i]);
0028 }
0029
0030 CSCSegment::~CSCSegment() {}
0031
0032 std::vector<const TrackingRecHit*> CSCSegment::recHits() const {
0033   std::vector<const TrackingRecHit*> pointersOfRecHits;
0034   for (std::vector<CSCRecHit2D>::const_iterator irh = theCSCRecHits.begin(); irh != theCSCRecHits.end(); ++irh) {
0035     pointersOfRecHits.push_back(&(*irh));
0036   }
0037   return pointersOfRecHits;
0038 }
0039
0040 std::vector<TrackingRecHit*> CSCSegment::recHits() {
0041   std::vector<TrackingRecHit*> pointersOfRecHits;
0042   for (std::vector<CSCRecHit2D>::iterator irh = theCSCRecHits.begin(); irh != theCSCRecHits.end(); ++irh) {
0043     pointersOfRecHits.push_back(&(*irh));
0044   }
0045   return pointersOfRecHits;
0046 }
0047
0048 LocalError CSCSegment::localPositionError() const {
0049   return LocalError(theCovMatrix[2][2], theCovMatrix[2][3], theCovMatrix[3][3]);
0050 }
0051
0052 LocalError CSCSegment::localDirectionError() const {
0053   return LocalError(theCovMatrix[0][0], theCovMatrix[0][1], theCovMatrix[1][1]);
0054 }
0055
0056 AlgebraicVector CSCSegment::parameters() const {
0057   // For consistency with DT and what we require for the TrackingRecHit interface,
0058   // the order of the parameters in the returned vector should be (dx/dz, dy/dz, x, z)
0059
0060   AlgebraicVector result(4);
0061
0062   result[0] = theLocalDirection.x() / theLocalDirection.z();
0063   result[1] = theLocalDirection.y() / theLocalDirection.z();
0064   result[2] = theOrigin.x();
0065   result[3] = theOrigin.y();
0066
0067   return result;
0068 }
0069
0070 AlgebraicMatrix createStaticMatrix() {
0071   AlgebraicMatrix m(4, 5, 0);
0072   m[0][1] = 1;
0073   m[1][2] = 1;
0074   m[2][3] = 1;
0075   m[3][4] = 1;
0076   return m;
0077 }
0078
0079 static const AlgebraicMatrix theProjectionMatrix = createStaticMatrix();
0080
0081 AlgebraicMatrix CSCSegment::projectionMatrix() const { return theProjectionMatrix; }
0082
0083 void CSCSegment::setDuplicateSegments(std::vector<CSCSegment*>& duplicates) {
0084   theDuplicateSegments.clear();
0085   for (unsigned int i = 0; i < duplicates.size(); ++i) {
0086     theDuplicateSegments.push_back(*duplicates[i]);
0087     //avoid copying duplicates of duplicates of duplicates...
0088     theDuplicateSegments.back().theDuplicateSegments.resize(0);
0089   }
0090 }
0091
0092 bool CSCSegment::testSharesAllInSpecificRecHits(const std::vector<CSCRecHit2D>& specificRecHits_1,
0093                                                 const std::vector<CSCRecHit2D>& specificRecHits_2,
0094                                                 CSCRecHit2D::SharedInputType sharesInput) const {
0095   const std::vector<CSCRecHit2D>* rhContainer_1 = &specificRecHits_1;
0096   const std::vector<CSCRecHit2D>* rhContainer_2 = &specificRecHits_2;
0097   if (specificRecHits_1.size() > specificRecHits_2.size()) {
0098     rhContainer_2 = &specificRecHits_1;
0099     rhContainer_1 = &specificRecHits_2;
0100   }
0101   //
0102   bool shareConditionPassed = true;
0103   for (std::vector<CSCRecHit2D>::const_iterator itRH = rhContainer_1->begin(); itRH != rhContainer_1->end(); ++itRH) {
0104     const CSCRecHit2D* firstRecHit = &(*itRH);
0105     bool sharedHit = false;
0106     for (std::vector<CSCRecHit2D>::const_iterator itRH2 = rhContainer_2->begin(); itRH2 != rhContainer_2->end();
0107          ++itRH2) {
0108       if (itRH2->sharesInput(firstRecHit, sharesInput)) {
0109         sharedHit = true;
0110         break;
0111       }
0112     }
0113     if (!sharedHit) {
0114       shareConditionPassed = false;
0115       break;
0116     }
0117   }
0118   return shareConditionPassed;
0119 }
0120
0121 //bool CSCSegment::sharesRecHits(CSCSegment  & anotherSegment, CSCRecHit2D::SharedInputType sharesInput){
0122 // 2 tracks through a chamber leave 4 rechits per layer (2 strips x 2 wire groups)
0123 // this function finds segments sharing wires or strips (first the rechits by sharesInput() )
0124 // there could probably be more complicated cases with partial sharing (but this needs studies)
0125 //
0126 //return testSharesAllInSpecificRecHits( theCSCRecHits , anotherSegment.specificRecHits(), sharesInput);
0127 //}
0128
0129 bool CSCSegment::sharesRecHits(const CSCSegment& anotherSegment, CSCRecHit2D::SharedInputType sharesInput) const {
0130   return testSharesAllInSpecificRecHits(theCSCRecHits, anotherSegment.specificRecHits(), sharesInput);
0131 }
0132
0133 //
0134 bool CSCSegment::sharesRecHits(const CSCSegment& anotherSegment) const {
0135   if (testSharesAllInSpecificRecHits(theCSCRecHits, anotherSegment.specificRecHits(), CSCRecHit2D::someWires) &&
0136       testSharesAllInSpecificRecHits(theCSCRecHits, anotherSegment.specificRecHits(), CSCRecHit2D::someStrips)) {
0137     return true;
0138   } else {
0139     return false;
0140   }
0141 }
0142 //
0143
0144 float CSCSegment::time() const {
0145   float averageTime = 0;
0146   std::vector<float> wireTimes;
0147   for (std::vector<CSCRecHit2D>::const_iterator itRH = theCSCRecHits.begin(); itRH != theCSCRecHits.end(); ++itRH) {
0148     const CSCRecHit2D* recHit = &(*itRH);
0149     averageTime += recHit->tpeak();
0150     averageTime += recHit->wireTime();
0151     wireTimes.push_back(recHit->wireTime());
0152   }
0153   averageTime = averageTime / (2 * theCSCRecHits.size());
0154
0155   //The wire times have a long tail that has to be pruned.  The strip times (tpeak) are fine
0156   bool modified = true;
0157   while (modified) {
0158     modified = false;
0159     double maxDiff = -1;
0160     std::vector<float>::iterator maxHit;
0161     for (std::vector<float>::iterator itWT = wireTimes.begin(); itWT != wireTimes.end(); ++itWT) {
0162       float diff = fabs(*itWT - averageTime);
0163       if (diff > maxDiff) {
0164         maxDiff = diff;
0165         maxHit = itWT;
0166       }
0167     }
0168     if (maxDiff > 26) {
0169       int N = theCSCRecHits.size() + wireTimes.size();
0170       averageTime = (averageTime * N - (*maxHit)) / (N - 1);
0171       wireTimes.erase(maxHit);
0172       modified = true;
0173     }
0174   }
0175   return averageTime;
0176 }
0177
0178 //
0179 void CSCSegment::print() const { std::cout << *this << std::endl; }
0180
0181 std::ostream& operator<<(std::ostream& os, const CSCSegment& seg) {
0182   os << "CSCSegment: local pos = " << seg.localPosition() << " posErr = (" << sqrt(seg.localPositionError().xx()) << ","
0183      << sqrt(seg.localPositionError().yy()) << "0,)\n"
0184      << "            dir = " << seg.localDirection() << " dirErr = (" << sqrt(seg.localDirectionError().xx()) << ","
0185      << sqrt(seg.localDirectionError().yy()) << "0,)\n"
0186      << "            chi2/ndf = " << seg.chi2() / double(seg.degreesOfFreedom())
0187      << " #rechits = " << seg.specificRecHits().size() << " ME1/1a duplicates : " << seg.duplicateSegments().size();
0188   return os;
0189 }
0190
0191 /*
0192 const CSCChamber* CSCSegment::chamber() const { return theChamber; }
0193 */