CSCSegment/src/CSCSegAlgoPreClustering.cc

0001 /**
0002  * \file CSCSegAlgoPreClustering.cc
0003  *
0004  *  \authors: S. Stoynev  - NU
0005  *            I. Bloch    - FNAL
0006  *            E. James    - FNAL
0007  *
0008  * ported by  D. Fortin   - UCR
0009  *
0010  * See header file for description.
0011  */
0012
0013 #include "RecoLocalMuon/CSCSegment/src/CSCSegAlgoPreClustering.h"
0014
0015 #include "Geometry/CSCGeometry/interface/CSCLayer.h"
0016 #include <Geometry/CSCGeometry/interface/CSCChamber.h>
0017 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0018
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0021
0022 #include <algorithm>
0023 #include <cmath>
0024 #include <iostream>
0025 #include <string>
0026
0027 /* Constructor
0028  *
0029  */
0030 CSCSegAlgoPreClustering::CSCSegAlgoPreClustering(const edm::ParameterSet& ps) {
0031   dXclusBoxMax = ps.getParameter<double>("dXclusBoxMax");
0032   dYclusBoxMax = ps.getParameter<double>("dYclusBoxMax");
0033   debug = ps.getUntrackedParameter<bool>("CSCSegmentDebug");
0034 }
0035
0036 /* Destructor:
0037  *
0038  */
0039 CSCSegAlgoPreClustering::~CSCSegAlgoPreClustering() {}
0040
0041 /* clusterHits
0042  *
0043  */
0044 std::vector<std::vector<const CSCRecHit2D*> > CSCSegAlgoPreClustering::clusterHits(const CSCChamber* aChamber,
0045                                                                                    const ChamberHitContainer& rechits) {
0046   theChamber = aChamber;
0047
0048   std::vector<ChamberHitContainer> rechits_clusters;  // this is a collection of groups of rechits
0049
0050   //float dXclus = 0.0;
0051   //float dYclus = 0.0;
0052   float dXclus_box = 0.0;
0053   float dYclus_box = 0.0;
0054
0055   std::vector<const CSCRecHit2D*> temp;
0056
0057   std::vector<ChamberHitContainer> seeds;
0058
0059   std::vector<float> running_meanX;
0060   std::vector<float> running_meanY;
0061
0062   std::vector<float> seed_minX;
0063   std::vector<float> seed_maxX;
0064   std::vector<float> seed_minY;
0065   std::vector<float> seed_maxY;
0066
0067   // split rechits into subvectors and return vector of vectors:
0068   // Loop over rechits
0069   // Create one seed per hit
0070   for (unsigned int i = 0; i < rechits.size(); ++i) {
0071     temp.clear();
0072
0073     temp.push_back(rechits[i]);
0074
0075     seeds.push_back(temp);
0076
0077     // First added hit in seed defines the mean to which the next hit is compared
0078     // for this seed.
0079
0080     running_meanX.push_back(rechits[i]->localPosition().x());
0081     running_meanY.push_back(rechits[i]->localPosition().y());
0082
0083     // set min/max X and Y for box containing the hits in the precluster:
0084     seed_minX.push_back(rechits[i]->localPosition().x());
0085     seed_maxX.push_back(rechits[i]->localPosition().x());
0086     seed_minY.push_back(rechits[i]->localPosition().y());
0087     seed_maxY.push_back(rechits[i]->localPosition().y());
0088   }
0089
0090   // merge clusters that are too close
0091   // measure distance between final "running mean"
0092   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
0093     for (size_t MMM = NNN + 1; MMM < seeds.size(); ++MMM) {
0094       if (running_meanX[MMM] == 999999. || running_meanX[NNN] == 999999.) {
0095         std::cout << "We should never see this line now!!!" << std::endl;
0096         continue;  //skip seeds that have been used
0097       }
0098
0099       // calculate cut criteria for simple running mean distance cut:
0100       //dXclus = fabs(running_meanX[NNN] - running_meanX[MMM]);
0101       //dYclus = fabs(running_meanY[NNN] - running_meanY[MMM]);
0102
0103       // calculate minmal distance between precluster boxes containing the hits:
0104       if (running_meanX[NNN] > running_meanX[MMM])
0105         dXclus_box = seed_minX[NNN] - seed_maxX[MMM];
0106       else
0107         dXclus_box = seed_minX[MMM] - seed_maxX[NNN];
0108       if (running_meanY[NNN] > running_meanY[MMM])
0109         dYclus_box = seed_minY[NNN] - seed_maxY[MMM];
0110       else
0111         dYclus_box = seed_minY[MMM] - seed_maxY[NNN];
0112
0113       if (dXclus_box < dXclusBoxMax && dYclus_box < dYclusBoxMax) {
0114         // merge clusters!
0115         // merge by adding seed NNN to seed MMM and erasing seed NNN
0116
0117         // calculate running mean for the merged seed:
0118         running_meanX[MMM] = (running_meanX[NNN] * seeds[NNN].size() + running_meanX[MMM] * seeds[MMM].size()) /
0119                              (seeds[NNN].size() + seeds[MMM].size());
0120         running_meanY[MMM] = (running_meanY[NNN] * seeds[NNN].size() + running_meanY[MMM] * seeds[MMM].size()) /
0121                              (seeds[NNN].size() + seeds[MMM].size());
0122
0123         // update min/max X and Y for box containing the hits in the merged cluster:
0124         if (seed_minX[NNN] <= seed_minX[MMM])
0125           seed_minX[MMM] = seed_minX[NNN];
0126         if (seed_maxX[NNN] > seed_maxX[MMM])
0127           seed_maxX[MMM] = seed_maxX[NNN];
0128         if (seed_minY[NNN] <= seed_minY[MMM])
0129           seed_minY[MMM] = seed_minY[NNN];
0130         if (seed_maxY[NNN] > seed_maxY[MMM])
0131           seed_maxY[MMM] = seed_maxY[NNN];
0132
0133         // add seed NNN to MMM (lower to larger number)
0134         seeds[MMM].insert(seeds[MMM].end(), seeds[NNN].begin(), seeds[NNN].end());
0135
0136         // mark seed NNN as used (at the moment just set running mean to 999999.)
0137         running_meanX[NNN] = 999999.;
0138         running_meanY[NNN] = 999999.;
0139         // we have merged a seed (NNN) to the highter seed (MMM) - need to contimue to
0140         // next seed (NNN+1)
0141         break;
0142       }
0143     }
0144   }
0145
0146   // hand over the final seeds to the output
0147   // would be more elegant if we could do the above step with
0148   // erasing the merged ones, rather than the
0149   for (size_t NNN = 0; NNN < seeds.size(); ++NNN) {
0150     if (running_meanX[NNN] == 999999.)
0151       continue;  //skip seeds that have been marked as used up in merging
0152     rechits_clusters.push_back(seeds[NNN]);
0153     mean_x = running_meanX[NNN];
0154     mean_y = running_meanY[NNN];
0155     err_x =
0156         (seed_maxX[NNN] - seed_minX[NNN]) / 3.464101615;  // use box size divided by sqrt(12) as position error estimate
0157     err_y =
0158         (seed_maxY[NNN] - seed_minY[NNN]) / 3.464101615;  // use box size divided by sqrt(12) as position error estimate
0159   }
0160
0161   return rechits_clusters;
0162 }