packages/MvtxPrototype2/MvtxPrototype2Clusterizer.cc

0001 #include "MvtxPrototype2Clusterizer.h"
0002 #include "MvtxPrototype2Geom.h"
0003
0004 //#include <g4detectors/PHG4CylinderCellGeom.h>
0005 //#include <g4detectors/PHG4CylinderCellGeomContainer.h>
0006 //#include <g4detectors/PHG4CylinderGeom.h>
0007 //#include <g4detectors/PHG4CylinderGeomContainer.h>
0008
0009 #include <mvtx/MvtxDefs.h>
0010 #include <mvtx/MvtxHit.h>
0011
0012 #include <trackbase/TrkrClusterContainer.h>
0013 #include <trackbase/TrkrClusterv1.h>
0014 #include <trackbase/TrkrHitSet.h>
0015 #include <trackbase/TrkrHitSetContainer.h>
0016 #include <trackbase/TrkrClusterHitAssoc.h>
0017
0018 #include <fun4all/Fun4AllReturnCodes.h>
0019 #include <phool/PHCompositeNode.h>
0020 #include <phool/PHIODataNode.h>
0021 #include <phool/PHNodeIterator.h>
0022 #include <phool/getClass.h>
0023
0024 #include <boost/format.hpp>
0025 #include <boost/tuple/tuple.hpp>
0026
0027 #include <TMatrixF.h>
0028 #include <TVector3.h>
0029
0030 #define BOOST_NO_HASH  // Our version of boost.graph is incompatible with GCC-4.3 w/o this flag
0031 #include <boost/bind.hpp>
0032 #include <boost/graph/adjacency_list.hpp>
0033 #include <boost/graph/connected_components.hpp>
0034 using namespace boost;
0035
0036 #include <cmath>
0037 #include <iostream>
0038 #include <stdexcept>
0039
0040 using namespace std;
0041
0042 static const float twopi = 2.0 * M_PI;
0043
0044 bool MvtxPrototype2Clusterizer::are_adjacent(const std::pair<TrkrDefs::hitkey, TrkrHit*> &lhs, const std::pair<TrkrDefs::hitkey, TrkrHit*> &rhs)
0045 {
0046   if (GetZClustering())
0047   {
0048     // column is first, row is second
0049     if (fabs( MvtxDefs::getCol(lhs.first) - MvtxDefs::getCol(rhs.first) ) <= 1)
0050     {
0051       if (fabs( MvtxDefs::getRow(lhs.first) - MvtxDefs::getRow(rhs.first) ) <= 1)
0052       {
0053         return true;
0054       }
0055     }
0056   }
0057   else
0058   {
0059     if (fabs( MvtxDefs::getCol(lhs.first) - MvtxDefs::getCol(rhs.first) ) == 0)
0060     {
0061       if (fabs( MvtxDefs::getRow(lhs.first) - MvtxDefs::getRow(rhs.first) ) <= 1)
0062       {
0063         return true;
0064       }
0065     }
0066   }
0067
0068   return false;
0069 }
0070
0071 MvtxPrototype2Clusterizer::MvtxPrototype2Clusterizer(const string &name)
0072   : SubsysReco(name)
0073   , m_hits(nullptr)
0074   , m_clusterlist(nullptr)
0075   , m_clusterhitassoc(nullptr)
0076   , m_makeZClustering(true)
0077 {
0078 }
0079
0080 int MvtxPrototype2Clusterizer::InitRun(PHCompositeNode *topNode)
0081 {
0082   //-----------------
0083   // Add Cluster Node
0084   //-----------------
0085
0086   PHNodeIterator iter(topNode);
0087
0088   // Looking for the DST node
0089   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0090   if (!dstNode)
0091   {
0092     cout << PHWHERE << "DST Node missing, doing nothing." << endl;
0093     return Fun4AllReturnCodes::ABORTRUN;
0094   }
0095   PHNodeIterator iter_dst(dstNode);
0096
0097   // Create the SVX node if required
0098   PHCompositeNode *svxNode = dynamic_cast<PHCompositeNode *>(iter_dst.findFirst("PHCompositeNode", "TRKR"));
0099   if (!svxNode)
0100   {
0101     svxNode = new PHCompositeNode("TRKR");
0102     dstNode->addNode(svxNode);
0103   }
0104
0105   // Create the Cluster node if required
0106   TrkrClusterContainer *trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
0107   if (!trkrclusters)
0108   {
0109     PHNodeIterator dstiter(dstNode);
0110     PHCompositeNode *DetNode =
0111       dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0112     if (!DetNode)
0113       {
0114     DetNode = new PHCompositeNode("TRKR");
0115     dstNode->addNode(DetNode);
0116       }
0117
0118     trkrclusters = new TrkrClusterContainer();
0119     PHIODataNode<PHObject> *TrkrClusterContainerNode =
0120       new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
0121     DetNode->addNode(TrkrClusterContainerNode);
0122   }
0123
0124   TrkrClusterHitAssoc *clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,"TRKR_CLUSTERHITASSOC");
0125   if(!clusterhitassoc)
0126     {
0127       PHNodeIterator dstiter(dstNode);
0128       PHCompositeNode *DetNode =
0129         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0130       if (!DetNode)
0131     {
0132       DetNode = new PHCompositeNode("TRKR");
0133       dstNode->addNode(DetNode);
0134     }
0135
0136       clusterhitassoc = new TrkrClusterHitAssoc();
0137       PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(clusterhitassoc, "TRKR_CLUSTERHITASSOC", "PHObject");
0138       DetNode->addNode(newNode);
0139     }
0140
0141   //----------------
0142   // Report Settings
0143   //----------------
0144
0145   if (Verbosity() > 0)
0146   {
0147     cout << "====================== MvtxPrototype2Clusterizer::InitRun() =====================" << endl;
0148     cout << " Z-dimension Clustering = " << boolalpha << m_makeZClustering << noboolalpha << endl;
0149     cout << "===========================================================================" << endl;
0150   }
0151
0152   return Fun4AllReturnCodes::EVENT_OK;
0153 }
0154
0155 int MvtxPrototype2Clusterizer::process_event(PHCompositeNode *topNode)
0156 {
0157   // get node containing the digitized hits
0158   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0159   if (!m_hits)
0160   {
0161     cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << endl;
0162     return Fun4AllReturnCodes::ABORTRUN;
0163   }
0164
0165   // get node for clusters
0166   m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0167   if (!m_clusterlist)
0168   {
0169     cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << endl;
0170     return Fun4AllReturnCodes::ABORTRUN;
0171   }
0172
0173   // get node for cluster hit associations
0174   m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0175   if (!m_clusterhitassoc)
0176   {
0177     cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << endl;
0178     return Fun4AllReturnCodes::ABORTRUN;
0179   }
0180
0181   // run clustering
0182   ClusterMvtx(topNode);
0183   PrintClusters(topNode);
0184
0185   // done
0186   return Fun4AllReturnCodes::EVENT_OK;
0187 }
0188
0189 void MvtxPrototype2Clusterizer::ClusterMvtx(PHCompositeNode *topNode)
0190 {
0191
0192     if (Verbosity() > 0)
0193         cout << "Entering MvtxPrototype2Clusterizer::ClusterMvtx " << endl;
0194
0195     //-----------
0196     // Clustering
0197     //-----------
0198
0199     // loop over each MvtxHitSet object (chip)
0200     TrkrHitSetContainer::ConstRange hitsetrange =
0201         m_hits->getHitSets(TrkrDefs::TrkrId::mvtxId);
0202     for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0203             hitsetitr != hitsetrange.second;
0204             ++hitsetitr)
0205     {
0206
0207         TrkrHitSet *hitset = hitsetitr->second;
0208
0209         if(Verbosity() > 10) cout << "MvtxPrototype2Clusterizer found hitsetkey " << hitsetitr->first << endl;
0210
0211         if (Verbosity() > 10)
0212             hitset->identify();
0213
0214         // fill a vector of hits to make things easier
0215         std::vector <std::pair< TrkrDefs::hitkey, TrkrHit*> > hitvec;
0216
0217         TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0218         for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0219                 hitr != hitrangei.second;
0220                 ++hitr)
0221         {
0222             hitvec.push_back(make_pair(hitr->first, hitr->second));
0223         }
0224         if (Verbosity() > 10)
0225             cout << "hitvec.size(): " << hitvec.size() << endl;
0226
0227         // do the clustering
0228         typedef adjacency_list<vecS, vecS, undirectedS> Graph;
0229         Graph G;
0230
0231         // loop over hits in this chip
0232         for (unsigned int i = 0; i < hitvec.size(); i++)
0233         {
0234             for (unsigned int j = 0; j < hitvec.size(); j++)
0235             {
0236                 if (are_adjacent(hitvec[i], hitvec[j]))
0237                     add_edge(i, j, G);
0238             }
0239         }
0240
0241         // Find the connections between the vertices of the graph (vertices are the rawhits,
0242         // connections are made when they are adjacent to one another)
0243         vector<int> component(num_vertices(G));
0244
0245         // this is the actual clustering, performed by boost
0246         connected_components(G, &component[0]);
0247
0248         // Loop over the components(hits) compiling a list of the
0249         // unique connected groups (ie. clusters).
0250         set<int> cluster_ids;  // unique components
0251         //multimap<int, pixel> clusters;
0252         multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*> >  clusters;
0253         for (unsigned int i = 0; i < component.size(); i++)
0254         {
0255             cluster_ids.insert(component[i]);
0256             clusters.insert(make_pair(component[i], hitvec[i]));
0257         }
0258
0259         // loop over the componenets and make clusters
0260         for (set<int>::iterator clusiter = cluster_ids.begin(); clusiter != cluster_ids.end(); ++clusiter)
0261         {
0262             int clusid = *clusiter;
0263             pair<multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*>>::iterator,
0264                 multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*>>::iterator>  clusrange = clusters.equal_range(clusid);
0265             multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*>>::iterator mapiter = clusrange.first;
0266
0267             if (Verbosity() > 2)
0268                 cout << "Filling cluster id " << clusid << endl;
0269
0270             // make the cluster directly in the node tree
0271             TrkrDefs::cluskey ckey = MvtxDefs::genClusKey(hitset->getHitSetKey(), clusid);
0272             TrkrClusterv1 *clus = static_cast<TrkrClusterv1 *>((m_clusterlist->findOrAddCluster(ckey))->second);
0273
0274             // we need the geometry object for this layer to get the global positions
0275             //int layer = TrkrDefs::getLayer(ckey);
0276             //int stave =  MvtxDefs::getStaveId(ckey);
0277             //int chip = MvtxDefs::getChipId(ckey);
0278
0279             // determine the size of the cluster in phi and z
0280             set<int> phibins;
0281             set<int> zbins;
0282
0283             // determine the cluster position...
0284             // need to define x, y, z coordinate
0285             double xsum = 0.0;
0286             double ysum = 0.0;
0287             double zsum = 0.0;
0288             unsigned nhits = 0;
0289
0290             double clusx = NAN;
0291             double clusy = NAN;
0292             double clusz = NAN;
0293
0294             for (mapiter = clusrange.first; mapiter != clusrange.second; ++mapiter)
0295             {
0296                 // size
0297                 int col =  MvtxDefs::getCol( (mapiter->second).first);
0298                 int row = MvtxDefs::getRow( (mapiter->second).first);
0299                 zbins.insert(col);
0300                 phibins.insert(row);
0301         double gloCoord[3];
0302         (MvtxPrototype2Geom::Instance())->detectorToGlobal(hitset->getHitSetKey(), (mapiter->second).first, gloCoord);
0303
0304                 // zsum += col;
0305                 // xsum += row;
0306                 // ysum += 0;
0307         xsum += gloCoord[0];
0308         ysum += gloCoord[1];
0309         zsum += gloCoord[2];
0310
0311                 // add the association between this cluster key and this hitkey to the table
0312                 m_clusterhitassoc->addAssoc(ckey, mapiter->second.first);
0313
0314                 ++nhits;
0315             }//mapiter
0316
0317             //float thickness = 50.e-4/28e-4; //sensor thickness converted to pixel size
0318       float thickness = SegmentationAlpide::SensorLayerThicknessEff;
0319             float phisize = phibins.size();
0320             float zsize = zbins.size();
0321
0322             // This is the local position
0323             // clusx = xsum/nhits + 0.5;
0324             // clusy = ysum/nhits + 0.5;
0325             // clusz = zsum/nhits + 0.5;
0326       clusx = xsum/nhits;
0327       clusy = ysum/nhits;
0328       clusz = zsum/nhits;
0329             clus->setAdc(nhits);
0330
0331             /*
0332             cout << "new mvtx clusterizer layer: " << layer << " stave: " << stave << " chip: " << chip
0333                 << " x: " << clusx << " y: " << clusy << " z: " << clusz << endl;
0334             */
0335
0336             clus->setPosition(0, clusx);
0337             clus->setPosition(1, clusy);
0338             clus->setPosition(2, clusz);
0339             //clus->setLocal();
0340       clus->setGlobal();
0341
0342             float invsqrt12 = 1.0 / sqrt(12.0);
0343
0344       phisize *= SegmentationAlpide::PitchRow;
0345       zsize   *= SegmentationAlpide::PitchCol;
0346             TMatrixF DIM(3, 3);
0347             DIM[0][0] = pow(0.5 * phisize, 2);
0348             DIM[0][1] = 0.0;
0349             DIM[0][2] = 0.0;
0350             DIM[1][0] = 0.0;
0351             DIM[1][1] = pow(0.5 * thickness, 2);
0352             DIM[1][2] = 0.0;
0353             DIM[2][0] = 0.0;
0354             DIM[2][1] = 0.0;
0355             DIM[2][2] = pow(0.5 * zsize, 2);
0356
0357             TMatrixF ERR(3, 3);
0358             ERR[0][0] = pow(0.5 * phisize * invsqrt12, 2);
0359             ERR[0][1] = 0.0;
0360             ERR[0][2] = 0.0;
0361             ERR[1][0] = 0.0;
0362             ERR[1][1] = pow(0.5 * thickness * invsqrt12, 2);
0363             ERR[1][2] = 0.0;
0364             ERR[2][0] = 0.0;
0365             ERR[2][1] = 0.0;
0366             ERR[2][2] = pow(0.5 * zsize * invsqrt12, 2);
0367
0368             clus->setSize( 0 , 0 , DIM[0][0] );
0369             clus->setSize( 0 , 1 , DIM[0][1] );
0370             clus->setSize( 0 , 2 , DIM[0][2] );
0371             clus->setSize( 1 , 0 , DIM[1][0] );
0372             clus->setSize( 1 , 1 , DIM[1][1] );
0373             clus->setSize( 1 , 2 , DIM[1][2] );
0374             clus->setSize( 2 , 0 , DIM[2][0] );
0375             clus->setSize( 2 , 1 , DIM[2][1] );
0376             clus->setSize( 2 , 2 , DIM[2][2] );
0377
0378             clus->setError( 0 , 0 , ERR[0][0] );
0379             clus->setError( 0 , 1 , ERR[0][1] );
0380             clus->setError( 0 , 2 , ERR[0][2] );
0381             clus->setError( 1 , 0 , ERR[1][0] );
0382             clus->setError( 1 , 1 , ERR[1][1] );
0383             clus->setError( 1 , 2 , ERR[1][2] );
0384             clus->setError( 2 , 0 , ERR[2][0] );
0385             clus->setError( 2 , 1 , ERR[2][1] );
0386             clus->setError( 2 , 2 , ERR[2][2] );
0387
0388             if (Verbosity() > 2)
0389                 clus->identify();
0390
0391         }
0392     }
0393
0394     if(Verbosity() > 5)
0395     {
0396         // check that the associations were written correctly
0397         m_clusterhitassoc->identify();
0398     }
0399
0400   return;
0401 }
0402
0403 void MvtxPrototype2Clusterizer::PrintClusters(PHCompositeNode *topNode)
0404 {
0405   if (Verbosity() >= 1)
0406   {
0407     TrkrClusterContainer *clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0408     if (!clusterlist) return;
0409
0410     cout << "================= Aftyer MvtxPrototype2Clusterizer::process_event() ====================" << endl;
0411
0412     cout << " There are " << clusterlist->size() << " clusters recorded: " << endl;
0413
0414     clusterlist->identify();
0415
0416     cout << "===========================================================================" << endl;
0417   }
0418
0419   return;
0420 }