packages/MvtxPrototype1/MvtxClusterizer.C

0001 #include "MvtxClusterizer.h"
0002
0003 #include "mvtx/MvtxDefUtil.h"
0004 #include "mvtx/MvtxHitSetv1.h"
0005
0006 #include <trackbase/TrkrHitSetContainer.h>
0007 #include <trackbase/TrkrHitSetv1.h>
0008 #include <trackbase/TrkrClusterContainer.h>
0009 #include <trackbase/TrkrClusterv1.h>
0010
0011
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <phool/PHCompositeNode.h>
0014 #include <phool/PHIODataNode.h>
0015 #include <phool/PHNodeIterator.h>
0016 #include <phool/getClass.h>
0017 #include <g4detectors/PHG4CellContainer.h>
0018 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0019 #include <g4detectors/PHG4CylinderGeomContainer.h>
0020 #include <g4detectors/PHG4CylinderGeom.h>
0021 #include <g4detectors/PHG4CylinderGeom_MAPS.h>
0022 #include <g4detectors/PHG4CylinderGeom_Siladders.h>
0023 #include <g4detectors/PHG4Cell.h>
0024 #include <g4detectors/PHG4CylinderCellGeom.h>
0025
0026 #include <boost/tuple/tuple.hpp>
0027 #include <boost/format.hpp>
0028
0029 #include <TMatrixF.h>
0030 #include <TVector3.h>
0031
0032 #define BOOST_NO_HASH // Our version of boost.graph is incompatible with GCC-4.3 w/o this flag
0033 #include <boost/bind.hpp>
0034 #include <boost/graph/adjacency_list.hpp>
0035 #include <boost/graph/connected_components.hpp>
0036 using namespace boost;
0037
0038 #include <iostream>
0039 #include <stdexcept>
0040 #include <cmath>
0041
0042 using namespace std;
0043
0044 static const float twopi = 2.0 * M_PI;
0045
0046 bool MvtxClusterizer::are_adjacent(const pixel lhs,
0047                                    const pixel rhs )
0048 {
0049
0050   if ( GetZClustering() )
0051   {
0052     // column is first, row is second
0053     if ( fabs(lhs.first - rhs.first) <= 1 )
0054     {
0055       if ( fabs(lhs.second - rhs.second) <= 1 )
0056       {
0057         return true;
0058       }
0059     }
0060   }
0061   else
0062   {
0063     if ( fabs(lhs.first - rhs.first) == 0 )
0064     {
0065       if ( fabs(lhs.second - rhs.second) <= 1 )
0066       {
0067         return true;
0068       }
0069     }
0070   }
0071
0072   return false;
0073 }
0074
0075
0076 MvtxClusterizer::MvtxClusterizer(const string &name) :
0077   SubsysReco(name),
0078   hits_(NULL),
0079   clusterlist_(NULL),
0080   make_z_clustering_(true),
0081   _timer(PHTimeServer::get()->insert_new(name))
0082 {
0083
0084 }
0085
0086 int MvtxClusterizer::InitRun(PHCompositeNode* topNode)
0087 {
0088
0089   //-----------------
0090   // Add Cluster Node
0091   //-----------------
0092
0093   PHNodeIterator iter(topNode);
0094
0095   // Looking for the DST node
0096   PHCompositeNode *dstNode
0097     = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0098   if (!dstNode) {
0099     cout << PHWHERE << "DST Node missing, doing nothing." << endl;
0100     return Fun4AllReturnCodes::ABORTRUN;
0101   }
0102   PHNodeIterator iter_dst(dstNode);
0103
0104   // Create the SVX node if required
0105   PHCompositeNode* svxNode
0106     = dynamic_cast<PHCompositeNode*>(iter_dst.findFirst("PHCompositeNode", "TRKR"));
0107   if (!svxNode) {
0108     svxNode = new PHCompositeNode("TRKR");
0109     dstNode->addNode(svxNode);
0110   }
0111
0112   // Create the Cluster node if required
0113   TrkrClusterContainer *trkrclusters
0114     = findNode::getClass<TrkrClusterContainer>(dstNode, "TrkrClusterContainer");
0115   if (!trkrclusters) {
0116     trkrclusters = new TrkrClusterContainer();
0117     PHIODataNode<PHObject> *TrkrClusterContainerNode =
0118       new PHIODataNode<PHObject>(trkrclusters, "TrkrClusterContainer", "PHObject");
0119     svxNode->addNode(TrkrClusterContainerNode);
0120   }
0121
0122   //----------------
0123   // Report Settings
0124   //----------------
0125
0126   if (verbosity > 0)
0127   {
0128     cout << "====================== MvtxClusterizer::InitRun() =====================" << endl;
0129     cout << " Z-dimension Clustering = " << boolalpha << make_z_clustering_ << noboolalpha << endl;
0130     cout << "===========================================================================" << endl;
0131   }
0132
0133   return Fun4AllReturnCodes::EVENT_OK;
0134 }
0135
0136 int MvtxClusterizer::process_event(PHCompositeNode *topNode) {
0137
0138   _timer.get()->restart();
0139
0140   // get node containing the digitized hits
0141   hits_ = findNode::getClass<TrkrHitSetContainer>(topNode, "TrkrHitSetContainer");
0142   if (!hits_)
0143   {
0144     cout << PHWHERE << "ERROR: Can't find node TrkrHitSetContainer" << endl;
0145     return Fun4AllReturnCodes::ABORTRUN;
0146   }
0147
0148   // get node for clusters
0149   clusterlist_ = findNode::getClass<TrkrClusterContainer>(topNode, "TrkrClusterContainer");
0150   if (!clusterlist_)
0151   {
0152     cout << PHWHERE << " ERROR: Can't find TrkrClusterContainer." << endl;
0153     return Fun4AllReturnCodes::ABORTRUN;
0154   }
0155     clusterlist_->Reset();
0156
0157   // run clustering
0158     ClusterMvtx(topNode);
0159   PrintClusters(topNode);
0160
0161   // done
0162   _timer.get()->stop();
0163   return Fun4AllReturnCodes::EVENT_OK;
0164 }
0165
0166 void MvtxClusterizer::ClusterMvtx(PHCompositeNode *topNode) {
0167
0168   if (verbosity > 0)
0169     cout << "Entering MvtxClusterizer::ClusterMvtx " << endl;
0170
0171   //----------
0172   // Get Nodes
0173   //----------
0174
0175   //-----------
0176   // Clustering
0177   //-----------
0178
0179   // loop over each MvtxHitSet object (chip)
0180   TrkrHitSetContainer::ConstRange hitsetrange =
0181     //hits_->GetHitSets(TrkrDefs::TRKRID::mvtx_id);
0182     hits_->GetHitSets();
0183   for ( TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0184         hitsetitr != hitsetrange.second;
0185         ++hitsetitr )
0186   {
0187
0188     MvtxHitSetv1* hitset = static_cast<MvtxHitSetv1*>(hitsetitr->second);
0189         if (verbosity>2)
0190             hitset->identify();
0191
0192     TrkrDefUtil trkrutil;
0193     MvtxDefUtil mvtxutil;
0194
0195
0196     // fill a vector of hits to make things easier
0197     // D. McGlinchey - this probably isn't necessary. Use iterator directly?
0198     std::vector<pixel> hitvec;
0199     MvtxHitSetv1::ConstRange hitrangei = hitset->GetHits();
0200     for ( MvtxHitSetv1::ConstIterator hitr = hitrangei.first;
0201           hitr != hitrangei.second;
0202           ++hitr)
0203     {
0204       hitvec.push_back(make_pair(hitr->first, hitr->second));
0205     }
0206         if ( verbosity>2 )
0207             cout << "hitvec.size(): " << hitvec.size() << endl;
0208
0209     // do the clustering
0210     typedef adjacency_list <vecS, vecS, undirectedS> Graph;
0211     Graph G;
0212
0213     // loop over hits in this chip
0214     for ( unsigned int i = 0; i < hitvec.size(); i++)
0215     {
0216       for ( unsigned int j = 0; j < hitvec.size(); j++)
0217       {
0218         if ( are_adjacent(hitvec[i], hitvec[j]) )
0219           add_edge(i, j, G);
0220       }
0221     }
0222
0223
0224     // Find the connections between the vertices of the graph (vertices are the rawhits,
0225     // connections are made when they are adjacent to one another)
0226     vector<int> component(num_vertices(G));
0227
0228     // this is the actual clustering, performed by boost
0229     connected_components(G, &component[0]);
0230
0231     // Loop over the components(hit cells) compiling a list of the
0232     // unique connected groups (ie. clusters).
0233     set<int> cluster_ids; // unique components
0234     multimap<int, pixel> clusters;
0235     for (unsigned int i = 0; i < component.size(); i++)
0236     {
0237       cluster_ids.insert( component[i] );
0238       clusters.insert( make_pair(component[i], hitvec[i]) );
0239     }
0240
0241     // loop over the componenets and make clusters
0242     for (set<int>::iterator clusiter = cluster_ids.begin();
0243          clusiter != cluster_ids.end();
0244          clusiter++ )
0245     {
0246
0247       int clusid = *clusiter;
0248       pair< multimap<int, pixel>::iterator,
0249             multimap<int, pixel>::iterator> clusrange = clusters.equal_range(clusid);
0250
0251       multimap<int, pixel>::iterator mapiter = clusrange.first;
0252
0253       if (verbosity > 2)
0254         cout << "Filling cluster id " << clusid << endl;
0255
0256       // make cluster
0257       TrkrDefs::cluskey ckey = mvtxutil.GenClusKey(hitset->GetHitSetKey(),clusid);
0258       TrkrClusterv1* clus = static_cast<TrkrClusterv1*>((clusterlist_->FindOrAddCluster(ckey))->second);
0259
0260       // determine the size of the cluster in phi and z
0261       set<int> phibins;
0262       set<int> zbins;
0263
0264       // determine the cluster position...
0265       double xsum = 0.0;
0266       double ysum = 0.0;
0267       double zsum = 0.0;
0268       unsigned nhits = 0;
0269
0270
0271       for (mapiter = clusrange.first; mapiter != clusrange.second; mapiter++ )
0272       {
0273
0274         // size
0275         zbins.insert((mapiter->second).first);
0276         phibins.insert((mapiter->second).second);
0277
0278         // find the center of the pixel in local coords
0279         xsum += (mapiter->second).second;
0280         zsum += (mapiter->second).first + 0.5;
0281
0282         ++nhits;
0283       } //mapitr
0284
0285       float thickness = 50.e-4/28e-4;
0286       float phisize = phibins.size();
0287       float zsize = zbins.size();
0288
0289       double clusx = NAN;
0290       double clusy = NAN;
0291       double clusz = NAN;
0292
0293       clusx = xsum / nhits;
0294       clusy = ysum / nhits;
0295       clusz = zsum / nhits;
0296
0297       clus->SetPosition(0, clusx);
0298       clus->SetPosition(1, clusy);
0299       clus->SetPosition(2, clusz);
0300       clus->SetLocal();
0301
0302       clus->SetAdc(nhits);
0303
0304       float invsqrt12 = 1.0 / sqrt(12.0);
0305
0306       TMatrixF DIM(3, 3);
0307       DIM[0][0] = pow(0.5 * phisize, 2);
0308
0309       DIM[0][1] = 0.0;
0310       DIM[0][2] = 0.0;
0311       DIM[1][0] = 0.0;
0312       DIM[1][1] = pow(0.5 * thickness, 2);
0313       DIM[1][2] = 0.0;
0314       DIM[2][0] = 0.0;
0315       DIM[2][1] = 0.0;
0316       DIM[2][2] = pow(0.5 * zsize, 2);
0317
0318       TMatrixF ERR(3, 3);
0319       ERR[0][0] = pow(0.5 * phisize * invsqrt12, 2);
0320       ERR[0][1] = 0.0;
0321       ERR[0][2] = 0.0;
0322       ERR[1][0] = 0.0;
0323       ERR[1][1] = pow(0.5 * thickness * invsqrt12, 2);
0324       ERR[1][2] = 0.0;
0325       ERR[2][0] = 0.0;
0326       ERR[2][1] = 0.0;
0327       ERR[2][2] = pow(0.5 * zsize * invsqrt12, 2);
0328
0329
0330       clus->SetSize( 0 , 0 , DIM[0][0] );
0331       clus->SetSize( 0 , 1 , DIM[0][1] );
0332       clus->SetSize( 0 , 2 , DIM[0][2] );
0333       clus->SetSize( 1 , 0 , DIM[1][0] );
0334       clus->SetSize( 1 , 1 , DIM[1][1] );
0335       clus->SetSize( 1 , 2 , DIM[1][2] );
0336       clus->SetSize( 2 , 0 , DIM[2][0] );
0337       clus->SetSize( 2 , 1 , DIM[2][1] );
0338       clus->SetSize( 2 , 2 , DIM[2][2] );
0339
0340       clus->SetError( 0 , 0 , ERR[0][0] );
0341       clus->SetError( 0 , 1 , ERR[0][1] );
0342       clus->SetError( 0 , 2 , ERR[0][2] );
0343       clus->SetError( 1 , 0 , ERR[1][0] );
0344       clus->SetError( 1 , 1 , ERR[1][1] );
0345       clus->SetError( 1 , 2 , ERR[1][2] );
0346       clus->SetError( 2 , 0 , ERR[2][0] );
0347       clus->SetError( 2 , 1 , ERR[2][1] );
0348       clus->SetError( 2 , 2 , ERR[2][2] );
0349
0350
0351       if ( verbosity > 2 )
0352         clus->identify();
0353
0354     } // clusitr
0355   } // hitsetitr
0356
0357   return;
0358 }
0359
0360
0361 void MvtxClusterizer::PrintClusters(PHCompositeNode * topNode)
0362 {
0363
0364   if (verbosity >= 1) {
0365
0366     TrkrClusterContainer *clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TrkrClusterContainer");
0367     if (!clusterlist) return;
0368
0369     cout << "================= MvtxClusterizer::process_event() ====================" << endl;
0370
0371     cout << " Found and recorded the following " << clusterlist->size() << " clusters: " << endl;
0372
0373     clusterlist->identify();
0374
0375     cout << "===========================================================================" << endl;
0376   }
0377
0378   return;
0379 }