ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MvtxClusterizer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MvtxClusterizer.cc
1 
7 #include "MvtxClusterizer.h"
8 #include "MvtxDefs.h"
9 #include "CylinderGeom_Mvtx.h"
10 
13 
16 #include <trackbase/TrkrDefs.h> // for hitkey, getLayer
17 #include <trackbase/TrkrHitv2.h>
18 #include <trackbase/TrkrHitSet.h>
21 
23 #include <fun4all/SubsysReco.h> // for SubsysReco
24 
25 #include <phool/PHCompositeNode.h>
26 #include <phool/PHIODataNode.h>
27 #include <phool/PHNode.h> // for PHNode
28 #include <phool/PHNodeIterator.h>
29 #include <phool/PHObject.h> // for PHObject
30 #include <phool/getClass.h>
31 #include <phool/phool.h> // for PHWHERE
32 
33 #include <TMatrixFfwd.h> // for TMatrixF
34 #include <TMatrixT.h> // for TMatrixT, operator*
35 #include <TMatrixTUtils.h> // for TMatrixTRow
36 #include <TVector3.h>
37 
38 #include <boost/graph/adjacency_list.hpp>
39 #include <boost/graph/connected_components.hpp>
40 
41 #include <array>
42 #include <cmath>
43 #include <cstdlib> // for exit
44 #include <iostream>
45 #include <map> // for multimap<>::iterator
46 #include <set> // for set, set<>::iterator
47 #include <string>
48 #include <vector> // for vector
49 
50 using namespace boost;
51 using namespace std;
52 
53 namespace
54 {
55 
57  template<class T>
58  inline constexpr T square( const T& x ) { return x*x; }
59 }
60 
61 bool MvtxClusterizer::are_adjacent(const std::pair<TrkrDefs::hitkey, TrkrHit*> &lhs, const std::pair<TrkrDefs::hitkey, TrkrHit*> &rhs)
62 {
63  if (GetZClustering())
64  {
65  // column is first, row is second
66  if (fabs( MvtxDefs::getCol(lhs.first) - MvtxDefs::getCol(rhs.first) ) <= 1)
67  {
68  if (fabs( MvtxDefs::getRow(lhs.first) - MvtxDefs::getRow(rhs.first) ) <= 1)
69  {
70  return true;
71  }
72  }
73  }
74  else
75  {
76  if (fabs( MvtxDefs::getCol(lhs.first) - MvtxDefs::getCol(rhs.first) ) == 0)
77  {
78  if (fabs( MvtxDefs::getRow(lhs.first) - MvtxDefs::getRow(rhs.first) ) <= 1)
79  {
80  return true;
81  }
82  }
83  }
84 
85  return false;
86 }
87 
89  : SubsysReco(name)
90  , m_hits(nullptr)
91  , m_clusterlist(nullptr)
92  , m_clusterhitassoc(nullptr)
93  , m_makeZClustering(true)
94 {
95 }
96 
98 {
99  //-----------------
100  // Add Cluster Node
101  //-----------------
102 
103  PHNodeIterator iter(topNode);
104 
105  // Looking for the DST node
106  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
107  if (!dstNode)
108  {
109  cout << PHWHERE << "DST Node missing, doing nothing." << endl;
111  }
112  PHNodeIterator iter_dst(dstNode);
113 
114  // Create the SVX node if required
115  PHCompositeNode *svxNode = dynamic_cast<PHCompositeNode *>(iter_dst.findFirst("PHCompositeNode", "TRKR"));
116  if (!svxNode)
117  {
118  svxNode = new PHCompositeNode("TRKR");
119  dstNode->addNode(svxNode);
120  }
121 
122  // Create the Cluster node if required
123  auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
124  if (!trkrclusters)
125  {
126  PHNodeIterator dstiter(dstNode);
127  PHCompositeNode *DetNode =
128  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
129  if (!DetNode)
130  {
131  DetNode = new PHCompositeNode("TRKR");
132  dstNode->addNode(DetNode);
133  }
134 
135  trkrclusters = new TrkrClusterContainerv3;
136  PHIODataNode<PHObject> *TrkrClusterContainerNode =
137  new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
138  DetNode->addNode(TrkrClusterContainerNode);
139  }
140 
141  auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,"TRKR_CLUSTERHITASSOC");
142  if(!clusterhitassoc)
143  {
144  PHNodeIterator dstiter(dstNode);
145  PHCompositeNode *DetNode =
146  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
147  if (!DetNode)
148  {
149  DetNode = new PHCompositeNode("TRKR");
150  dstNode->addNode(DetNode);
151  }
152 
153  clusterhitassoc = new TrkrClusterHitAssocv3;
154  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(clusterhitassoc, "TRKR_CLUSTERHITASSOC", "PHObject");
155  DetNode->addNode(newNode);
156  }
157 
158 
159  //----------------
160  // Report Settings
161  //----------------
162 
163  if (Verbosity() > 0)
164  {
165  cout << "====================== MvtxClusterizer::InitRun() =====================" << endl;
166  cout << " Z-dimension Clustering = " << boolalpha << m_makeZClustering << noboolalpha << endl;
167  cout << "===========================================================================" << endl;
168  }
169 
171 }
172 
174 {
175  // get node containing the digitized hits
176  m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
177  if (!m_hits)
178  {
179  cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << endl;
181  }
182 
183  // get node for clusters
184  m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
185  if (!m_clusterlist)
186  {
187  cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << endl;
189  }
190 
191  // get node for cluster hit associations
192  m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
193  if (!m_clusterhitassoc)
194  {
195  cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << endl;
197  }
198 
199  // run clustering
200  ClusterMvtx(topNode);
201  PrintClusters(topNode);
202 
203  // done
205 }
206 
208 {
209  if (Verbosity() > 0)
210  cout << "Entering MvtxClusterizer::ClusterMvtx " << endl;
211 
212  PHG4CylinderGeomContainer* geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
213  if (!geom_container) return;
214 
215  //-----------
216  // Clustering
217  //-----------
218 
219  // loop over each MvtxHitSet object (chip)
220  TrkrHitSetContainer::ConstRange hitsetrange =
222  for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
223  hitsetitr != hitsetrange.second;
224  ++hitsetitr)
225  {
226  TrkrHitSet *hitset = hitsetitr->second;
227 
228  if(Verbosity() > 1) cout << "MvtxClusterizer found hitsetkey " << hitsetitr->first << endl;
229 
230  if (Verbosity() > 2)
231  hitset->identify();
232 
233  // fill a vector of hits to make things easier
234  std::vector <std::pair< TrkrDefs::hitkey, TrkrHit*> > hitvec;
235 
236  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
237  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
238  hitr != hitrangei.second;
239  ++hitr)
240  {
241  hitvec.push_back(make_pair(hitr->first, hitr->second));
242  }
243  if (Verbosity() > 2) cout << "hitvec.size(): " << hitvec.size() << endl;
244 
245  // do the clustering
246  typedef adjacency_list<vecS, vecS, undirectedS> Graph;
247  Graph G;
248 
249  // loop over hits in this chip
250  for (unsigned int i = 0; i < hitvec.size(); i++)
251  {
252  for (unsigned int j = 0; j < hitvec.size(); j++)
253  {
254  if (are_adjacent(hitvec[i], hitvec[j]))
255  add_edge(i, j, G);
256  }
257  }
258 
259  // Find the connections between the vertices of the graph (vertices are the rawhits,
260  // connections are made when they are adjacent to one another)
261  vector<int> component(num_vertices(G));
262 
263  // this is the actual clustering, performed by boost
264  connected_components(G, &component[0]);
265 
266  // Loop over the components(hits) compiling a list of the
267  // unique connected groups (ie. clusters).
268  set<int> cluster_ids; // unique components
269  //multimap<int, pixel> clusters;
270  multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*> > clusters;
271  for (unsigned int i = 0; i < component.size(); i++)
272  {
273  cluster_ids.insert(component[i]);
274  clusters.insert(make_pair(component[i], hitvec[i]));
275  }
276  // cout << "found cluster #: "<< clusters.size()<< endl;
277  // loop over the componenets and make clusters
278  for (set<int>::iterator clusiter = cluster_ids.begin(); clusiter != cluster_ids.end(); ++clusiter)
279  {
280  int clusid = *clusiter;
281  auto clusrange = clusters.equal_range(clusid);
282 
283  if (Verbosity() > 2) cout << "Filling cluster id " << clusid << " of " << std::distance(cluster_ids.begin(),clusiter )<< endl;
284 
285  // make the cluster directly in the node tree
286  auto ckey = MvtxDefs::genClusKey(hitset->getHitSetKey(), clusid);
287 
288  auto clus = std::make_unique<TrkrClusterv3>();
289  clus->setClusKey(ckey);
290 
291  // determine the size of the cluster in phi and z
292  set<int> phibins;
293  set<int> zbins;
294 
295  // determine the cluster position...
296  double locxsum = 0.;
297  double loczsum = 0.;
298  const unsigned int nhits = std::distance( clusrange.first, clusrange.second );
299 
300  double locclusx = NAN;
301  double locclusz = NAN;
302 
303  // we need the geometry object for this layer to get the global positions
304  int layer = TrkrDefs::getLayer(ckey);
305  auto layergeom = dynamic_cast<CylinderGeom_Mvtx *>(geom_container->GetLayerGeom(layer));
306  if (!layergeom)
307  exit(1);
308 
309  for ( auto mapiter = clusrange.first; mapiter != clusrange.second; ++mapiter)
310  {
311  // size
312  int col = MvtxDefs::getCol( (mapiter->second).first);
313  int row = MvtxDefs::getRow( (mapiter->second).first);
314  zbins.insert(col);
315  phibins.insert(row);
316 
317  // get local coordinates, in stae reference frame, for hit
318  auto local_coords = layergeom->get_local_coords_from_pixel(row,col);
319 
320  /*
321  manually offset position along y (thickness of the sensor),
322  to account for effective hit position in the sensor, resulting from diffusion.
323  Effective position corresponds to 1um above the middle of the sensor
324  */
325  local_coords.SetY( 1e-4 );
326 
327  // update cluster position
328  locxsum += local_coords.X();
329  loczsum += local_coords.Z();
330  // add the association between this cluster key and this hitkey to the table
331  m_clusterhitassoc->addAssoc(ckey, mapiter->second.first);
332 
333  } //mapiter
334 
335  // This is the local position
336  locclusx = locxsum / nhits;
337  locclusz = loczsum / nhits;
338 
339  clus->setAdc(nhits);
340 
341  const double pitch = layergeom->get_pixel_x();
342  const double length = layergeom->get_pixel_z();
343  const double phisize = phibins.size() * pitch;
344  const double zsize = zbins.size() * length;
345 
346  static const double invsqrt12 = 1./std::sqrt(12);
347 
348  // scale factors (phi direction)
349  /*
350  they corresponds to clusters of size (2,2), (2,3), (3,2) and (3,3) in phi and z
351  other clusters, which are very few and pathological, get a scale factor of 1
352  These scale factors are applied to produce cluster pulls with width unity
353  */
354 
355  double phierror = pitch * invsqrt12;
356 
357  static constexpr std::array<double, 7> scalefactors_phi = {{ 0.36, 0.6,0.37,0.49,0.4,0.37,0.33 }};
358  if(phibins.size() == 1 && zbins.size() == 1) phierror*=scalefactors_phi[0];
359  else if(phibins.size() == 2 && zbins.size() == 1) phierror*=scalefactors_phi[1];
360  else if(phibins.size() == 1 && zbins.size() == 2) phierror*=scalefactors_phi[2];
361  else if( phibins.size() == 2 && zbins.size() == 2 ) phierror*=scalefactors_phi[0];
362  else if( phibins.size() == 2 && zbins.size() == 3 ) phierror*=scalefactors_phi[1];
363  else if( phibins.size() == 3 && zbins.size() == 2 ) phierror*=scalefactors_phi[2];
364  else if( phibins.size() == 3 && zbins.size() == 3 ) phierror*=scalefactors_phi[3];
365 
366 
367  // scale factors (z direction)
368  /*
369  they corresponds to clusters of size (2,2), (2,3), (3,2) and (3,3) in z and phi
370  other clusters, which are very few and pathological, get a scale factor of 1
371  */
372  static constexpr std::array<double, 4> scalefactors_z = {{ 0.47, 0.48, 0.71, 0.55 }};
373  double zerror = length*invsqrt12;
374  if( zbins.size() == 2 && phibins.size() == 2 ) zerror*=scalefactors_z[0];
375  else if( zbins.size() == 2 && phibins.size() == 3 ) zerror*=scalefactors_z[1];
376  else if( zbins.size() == 3 && phibins.size() == 2 ) zerror*=scalefactors_z[2];
377  else if( zbins.size() == 3 && phibins.size() == 3 ) zerror*=scalefactors_z[3];
378 
379  if(Verbosity() > 0)
380  cout << " MvtxClusterizer: layer " << layer << " rad " << layergeom->get_radius() << " phibins " << phibins.size() << " pitch " << pitch << " phisize " << phisize
381  << " zbins " << zbins.size() << " length " << length << " zsize " << zsize << endl;
382 
383  clus->setLocalX(locclusx);
384  clus->setLocalY(locclusz);
386  clus->setActsLocalError(0,0,square(phierror));
387  clus->setActsLocalError(0,1,0.);
388  clus->setActsLocalError(1,0,0.);
389  clus->setActsLocalError(1,1,square(zerror));
390 
393  clus->setSubSurfKey(0);
394 
395  if (Verbosity() > 2)
396  clus->identify();
397 
398  m_clusterlist->addCluster(clus.release());
399 
400  } // clusitr
401  } // hitsetitr
402 
403  if(Verbosity() > 1)
404  {
405  // check that the associations were written correctly
407  }
408 
409  return;
410 }
411 
413 {
414  if (Verbosity() >= 1)
415  {
416  TrkrClusterContainer *clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
417  if (!clusterlist) return;
418 
419  cout << "================= Aftyer MvtxClusterizer::process_event() ====================" << endl;
420 
421  cout << " There are " << clusterlist->size() << " clusters recorded: " << endl;
422 
423  clusterlist->identify();
424 
425  cout << "===========================================================================" << endl;
426  }
427 
428  return;
429 }
430