packages/CaloReco/PHMakeGroups.h

0001 #ifndef CALORECO_PHMAKEGROUPS_H
0002 #define CALORECO_PHMAKEGROUPS_H
0003 // Requirements:
0004 //
0005 // the class type Hit needs to provide:
0006 //
0007 //   an operator< that sorts by ix and then by iz
0008 //   a function is_adjacent() which returns true if the argumment hit is adjacent
0009 //   a function ...() that returns true if the argument Hit is far enough away
0010 //      from this one to allow breaking out of the inner loop early
0011 //
0012
0013 #include <boost/bind/bind.hpp>
0014 #pragma GCC diagnostic push
0015 #pragma GCC diagnostic ignored "-Wshadow"
0016 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
0017 #include <boost/graph/adjacency_list.hpp>
0018 #include <boost/graph/connected_components.hpp>
0019 #pragma GCC diagnostic pop
0020
0021 #include <map>
0022 #include <vector>
0023
0024 template <class Hit>
0025 int PHMakeGroups(std::vector<Hit>& hits,
0026                  std::multimap<int, Hit>& groups)
0027 {
0028   typedef boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS> Graph;
0029
0030   Graph G;
0031
0032   // Process the hits in ix-then-iz order
0033   std::sort(hits.begin(), hits.end());
0034
0035   // TODO: Since the list is sorted by channel number, it should
0036   // be possible to terminate the inner loop if we find the next hit
0037   // is more than one Z channel away (in which case the subsequent ones
0038   // will be adjacent neither in X nor Z).
0039   for (unsigned int i = 0; i < hits.size(); i++)
0040   {
0041     for (unsigned int j = i + 1; j < hits.size(); j++)
0042     {
0043       if (hits[i].is_adjacent(hits[j])) add_edge(i, j, G);
0044     }
0045     add_edge(i, i, G);
0046   }
0047
0048   // Find the connections between the vertices of the graph (vertices are the rawhits,
0049   // connections are made when they are adjacent to one another)
0050   std::vector<int> component(num_vertices(G));
0051   // connected_components(G, &component[0]);
0052   connected_components(G, &component[0]);
0053   // std::cout << "Found " << num << " groups of hits" << std::endl;
0054
0055   // Loop over the components(vertices) compiling a list of the unique
0056   // connections (ie clusters).
0057   std::set<int> comps;  // Number of unique components
0058   for (unsigned int i = 0; i < component.size(); i++)
0059   {
0060     comps.insert(component[i]);
0061     groups.insert(std::make_pair(component[i], hits[i]));
0062   }
0063
0064   //       for(std::set<int>::const_iterator id=comps.begin(); id!=comps.end(); id++)
0065   //    {
0066   //      std::multimap<int,SvxRawhitAdapter>::const_iterator curr, last;
0067   //      boost::tie(curr,last) = groups[iread].equal_range(*id);
0068   //      std::cout << "Group " << *id << " has " << groups[iread].count(*id) << " Rawhits:" << std::endl;
0069   //      for ( ; curr!=last; curr++)
0070   //        {
0071   //          SvxRawhitAdapter h = curr->second;
0072   //          h.hit->print();
0073   //        }
0074   //    }
0075
0076   return 0;
0077 }
0078
0079 #endif