ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
InttClusterizer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file InttClusterizer.cc
1 #include "InttClusterizer.h"
2 #include "CylinderGeomIntt.h"
3 #include "InttDefs.h"
4 
7 #include <trackbase/TrkrDefs.h>
8 #include <trackbase/TrkrHitSet.h>
9 #include <trackbase/TrkrHitv2.h>
12 
15 
17 #include <fun4all/SubsysReco.h>
18 
19 #include <phool/PHCompositeNode.h>
20 #include <phool/PHIODataNode.h>
21 #include <phool/PHNode.h>
22 #include <phool/PHNodeIterator.h>
23 #include <phool/PHObject.h> // for PHObject
24 #include <phool/getClass.h>
25 #include <phool/phool.h>
26 
27 #include <TMatrixFfwd.h> // for TMatrixF
28 #include <TMatrixT.h> // for TMatrixT, operator*
29 #include <TMatrixTUtils.h> // for TMatrixTRow
30 
31 #include <boost/graph/adjacency_list.hpp>
32 #include <boost/graph/connected_components.hpp>
33 
34 #include <array>
35 #include <cmath>
36 #include <iostream>
37 #include <set>
38 #include <vector> // for vector
39 
40 using namespace boost;
41 using namespace std;
42 
43 namespace
44 {
45 
47  template<class T>
48  inline constexpr T square( const T& x ) { return x*x; }
49 }
50 
51 bool InttClusterizer::ladder_are_adjacent( const std::pair<TrkrDefs::hitkey, TrkrHit*> &lhs, const std::pair<TrkrDefs::hitkey, TrkrHit*> &rhs, const int layer)
52 {
53  if (get_z_clustering(layer))
54  {
55  if (fabs( InttDefs::getCol(lhs.first) - InttDefs::getCol(rhs.first) ) <= 1)
56  {
57  if (fabs( InttDefs::getRow(lhs.first) - InttDefs::getRow(rhs.first) ) <= 1)
58  {
59  return true;
60  }
61  }
62  }
63  else
64  if (fabs( InttDefs::getCol(lhs.first) - InttDefs::getCol(rhs.first) ) == 0)
65  {
66  if (fabs( InttDefs::getRow(lhs.first) - InttDefs::getRow(rhs.first) ) <= 1)
67  {
68  return true;
69  }
70  }
71 
72  return false;
73 }
74 
76  unsigned int /*min_layer*/,
77  unsigned int /*max_layer*/)
78  : SubsysReco(name)
79  , m_hits(nullptr)
80  , m_clusterlist(nullptr)
81  , m_clusterhitassoc(nullptr)
82  , _fraction_of_mip(0.5)
83  , _thresholds_by_layer()
84  , _make_z_clustering()
85  , _make_e_weights()
86 {
87 }
88 
90 {
91  /*
92  // get node containing the digitized hits
93  _hits = findNode::getClass<SvtxHitMap>(topNode, "SvtxHitMap");
94  if (!_hits)
95  {
96  cout << PHWHERE << "ERROR: Can't find node SvtxHitMap" << endl;
97  return Fun4AllReturnCodes::ABORTRUN;
98  }
99  */
100 
101  //-----------------
102  // Add Cluster Node
103  //-----------------
104 
105  PHNodeIterator iter(topNode);
106 
107  // Looking for the DST node
108  PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
109  if (!dstNode)
110  {
111  cout << PHWHERE << "DST Node missing, doing nothing." << endl;
113  }
114  PHNodeIterator iter_dst(dstNode);
115 
116  // Create the Cluster and association nodes if required
117  auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
118  if (!trkrclusters)
119  {
120  PHNodeIterator dstiter(dstNode);
121  PHCompositeNode *DetNode =
122  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
123  if (!DetNode)
124  {
125  DetNode = new PHCompositeNode("TRKR");
126  dstNode->addNode(DetNode);
127  }
128 
129  trkrclusters = new TrkrClusterContainerv3;
130  PHIODataNode<PHObject> *TrkrClusterContainerNode =
131  new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
132  DetNode->addNode(TrkrClusterContainerNode);
133  }
134 
135  auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,"TRKR_CLUSTERHITASSOC");
136  if(!clusterhitassoc)
137  {
138  PHNodeIterator dstiter(dstNode);
139  PHCompositeNode *DetNode =
140  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
141  if (!DetNode)
142  {
143  DetNode = new PHCompositeNode("TRKR");
144  dstNode->addNode(DetNode);
145  }
146 
147  clusterhitassoc = new TrkrClusterHitAssocv3;
148  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(clusterhitassoc, "TRKR_CLUSTERHITASSOC", "PHObject");
149  DetNode->addNode(newNode);
150  }
151 
152 
153  //---------------------
154  // Calculate Thresholds
155  //---------------------
156 
157  CalculateLadderThresholds(topNode);
158 
159  //----------------
160  // Report Settings
161  //----------------
162 
163  if (Verbosity() > 0)
164  {
165  cout << "====================== InttClusterizer::InitRun() =====================" << endl;
166  cout << " Fraction of expected thickness MIP energy = " << _fraction_of_mip << endl;
167  for (std::map<int, float>::iterator iter = _thresholds_by_layer.begin();
168  iter != _thresholds_by_layer.end();
169  ++iter)
170  {
171  cout << " Cluster Threshold in Layer #" << iter->first << " = " << 1.0e6 * iter->second << " keV" << endl;
172  }
173  for (std::map<int, bool>::iterator iter = _make_z_clustering.begin();
174  iter != _make_z_clustering.end();
175  ++iter)
176  {
177  cout << " Z-dimension Clustering in Layer #" << iter->first << " = " << boolalpha << iter->second << noboolalpha << endl;
178  }
179  for (std::map<int, bool>::iterator iter = _make_e_weights.begin();
180  iter != _make_e_weights.end();
181  ++iter)
182  {
183  cout << " Energy weighting clusters in Layer #" << iter->first << " = " << boolalpha << iter->second << noboolalpha << endl;
184  }
185  cout << "===========================================================================" << endl;
186  }
187 
189 }
190 
192 {
193 
194  // get node containing the digitized hits
195  m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
196  if (!m_hits)
197  {
198  cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << endl;
200  }
201 
202  // get node for clusters
203  m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
204  if (!m_clusterlist)
205  {
206  cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << endl;
208  }
209 
210  // get node for cluster hit associations
211  m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
212  if (!m_clusterhitassoc)
213  {
214  cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << endl;
216  }
217 
218 
219 
220  ClusterLadderCells(topNode);
221  PrintClusters(topNode);
222 
224 }
225 
227 {
228  /*
229  PHG4CellContainer* cells = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_INTT");
230  if (!cells) return;
231  */
232 
233  PHG4CylinderGeomContainer* geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
234  if (!geom_container) return;
235 
236  PHG4CylinderGeomContainer::ConstRange layerrange = geom_container->get_begin_end();
237  for (PHG4CylinderGeomContainer::ConstIterator layeriter = layerrange.first;
238  layeriter != layerrange.second;
239  ++layeriter)
240  {
241  // index mapping
242  int layer = layeriter->second->get_layer();
243 
244  // cluster threshold
245  float thickness = (layeriter->second)->get_thickness();
246  float threshold = _fraction_of_mip * 0.003876 * thickness;
247  _thresholds_by_layer.insert(std::make_pair(layer, threshold));
248 
249  // fill in a default z_clustering value if not present
250  if (_make_z_clustering.find(layer) == _make_z_clustering.end())
251  {
252  _make_z_clustering.insert(std::make_pair(layer, false));
253  }
254 
255  if (_make_e_weights.find(layer) == _make_e_weights.end())
256  {
257  _make_e_weights.insert(std::make_pair(layer, false));
258  }
259  }
260 
261  return;
262 }
263 
265 {
266  if (Verbosity() > 0)
267  cout << "Entering InttClusterizer::ClusterLadderCells " << endl;
268 
269  //----------
270  // Get Nodes
271  //----------
272 
273  // get the geometry node
274  PHG4CylinderGeomContainer* geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
275  if (!geom_container) return;
276 
277  //-----------
278  // Clustering
279  //-----------
280 
281  // loop over the InttHitSet objects
282  TrkrHitSetContainer::ConstRange hitsetrange =
284  for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
285  hitsetitr != hitsetrange.second;
286  ++hitsetitr)
287  {
288  // Each hitset contains only hits that are clusterizable - i.e. belong to a single sensor
289  TrkrHitSet *hitset = hitsetitr->second;
290 
291  if(Verbosity() > 1) cout << "InttClusterizer found hitsetkey " << hitsetitr->first << endl;
292  if (Verbosity() > 2)
293  hitset->identify();
294 
295  // we have a single hitset, get the info that identifies the sensor
296  int layer = TrkrDefs::getLayer(hitsetitr->first);
297  int ladder_z_index = InttDefs::getLadderZId(hitsetitr->first);
298 
299  // we will need the geometry object for this layer to get the global position
300  CylinderGeomIntt* geom = dynamic_cast<CylinderGeomIntt*>(geom_container->GetLayerGeom(layer));
301  float pitch = geom->get_strip_y_spacing();
302  float length = geom->get_strip_z_spacing();
303 
304  // fill a vector of hits to make things easier - gets every hit in the hitset
305  std::vector <std::pair< TrkrDefs::hitkey, TrkrHit*> > hitvec;
306  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
307  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
308  hitr != hitrangei.second;
309  ++hitr)
310  {
311  hitvec.push_back(make_pair(hitr->first, hitr->second));
312  }
313  if (Verbosity() > 2)
314  cout << "hitvec.size(): " << hitvec.size() << endl;
315 
316  typedef adjacency_list<vecS, vecS, undirectedS> Graph;
317  Graph G;
318 
319  // Find adjacent strips
320  for (unsigned int i = 0; i < hitvec.size(); i++)
321  {
322  for (unsigned int j = i + 1; j < hitvec.size(); j++)
323  {
324  if (ladder_are_adjacent(hitvec[i], hitvec[j], layer))
325  {
326  add_edge(i, j, G);
327  }
328  }
329 
330  add_edge(i, i, G);
331  }
332 
333  // Find the connections between the vertices of the graph (vertices are the rawhits,
334  // connections are made when they are adjacent to one another)
335  vector<int> component(num_vertices(G));
336 
337  // this is the actual clustering, performed by boost
338  connected_components(G, &component[0]);
339 
340  // Loop over the components(hit cells) compiling a list of the
341  // unique connected groups (ie. clusters).
342  set<int> cluster_ids; // unique components
343 
344  multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*> > clusters;
345  for (unsigned int i = 0; i < component.size(); i++)
346  {
347  cluster_ids.insert(component[i]); // one entry per unique cluster id
348  clusters.insert(make_pair(component[i], hitvec[i])); // multiple entries per unique cluster id
349  }
350 
351  // loop over the cluster ID's and make the clusters from the connected hits
352  for (set<int>::iterator clusiter = cluster_ids.begin(); clusiter != cluster_ids.end(); ++clusiter)
353  {
354  int clusid = *clusiter;
355  //cout << " intt clustering: add cluster number " << clusid << endl;
356  // get all hits for this cluster ID only
357  pair<multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*>>::iterator,
358  multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*>>::iterator> clusrange = clusters.equal_range(clusid);
359  multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*>>::iterator mapiter = clusrange.first;
360 
361  // make the cluster directly in the node tree
362  TrkrDefs::cluskey ckey = InttDefs::genClusKey(hitset->getHitSetKey(), clusid);
363  auto clus = std::make_unique<TrkrClusterv3>();
364  clus->setClusKey(ckey);
365 
366  if (Verbosity() > 2)
367  cout << "Filling cluster with key " << ckey << endl;
368 
369  // determine the size of the cluster in phi and z, useful for track fitting the cluster
370  set<int> phibins;
371  set<int> zbins;
372 
373  // determine the cluster position...
374  double xlocalsum = 0.0;
375  double ylocalsum = 0.0;
376  double zlocalsum = 0.0;
377  unsigned int clus_adc = 0.0;
378  unsigned nhits = 0;
379 
380  for (mapiter = clusrange.first; mapiter != clusrange.second; ++mapiter)
381  {
382  // mapiter->second.first is the hit key
383  //cout << " adding hitkey " << mapiter->second.first << endl;
384  int col = InttDefs::getCol( (mapiter->second).first);
385  int row = InttDefs::getRow( (mapiter->second).first);
386  zbins.insert(col);
387  phibins.insert(row);
388 
389  // mapiter->second.second is the hit
390  unsigned int hit_adc = (mapiter->second).second->getAdc();
391 
392  // now get the positions from the geometry
393  double local_hit_location[3] = {0., 0., 0.};
394  geom->find_strip_center_localcoords(ladder_z_index,
395  row, col,
396  local_hit_location);
397 
398  if (_make_e_weights[layer])
399  {
400  xlocalsum += local_hit_location[0] * (double) hit_adc;
401  ylocalsum += local_hit_location[1] * (double) hit_adc;
402  zlocalsum += local_hit_location[2] * (double) hit_adc;
403  }
404  else
405  {
406  xlocalsum += local_hit_location[0];
407  ylocalsum += local_hit_location[1];
408  zlocalsum += local_hit_location[2];
409  }
410 
411  clus_adc += hit_adc;
412  ++nhits;
413 
414  // add this cluster-hit association to the association map of (clusterkey,hitkey)
415  m_clusterhitassoc->addAssoc(ckey, mapiter->second.first);
416 
417  if (Verbosity() > 2) cout << " nhits = " << nhits << endl;
418  if (Verbosity() > 2)
419  {
420  cout << " From geometry object: hit x " << local_hit_location[0] << " hit y " << local_hit_location[1] << " hit z " << local_hit_location[2] << endl;
421  cout << " nhits " << nhits << " clusx = " << xlocalsum / nhits << " clusy " << ylocalsum / nhits << " clusz " << zlocalsum / nhits << " hit_adc " << hit_adc << endl;
422 
423  }
424  }
425 
426  static const float invsqrt12 = 1./sqrt(12);
427 
428  // scale factors (phi direction)
429  /*
430  they corresponds to clusters of size 1 and 2 in phi
431  other clusters, which are very few and pathological, get a scale factor of 1
432  These scale factors are applied to produce cluster pulls with width unity
433  */
434 
435  float phierror = pitch * invsqrt12;
436 
437  static constexpr std::array<double, 3> scalefactors_phi = {{ 0.85, 0.4, 0.33 }};
438  if( phibins.size() == 1 && layer < 5) phierror*=scalefactors_phi[0];
439  else if( phibins.size() == 2 && layer < 5) phierror*=scalefactors_phi[1];
440  else if( phibins.size() == 2 && layer > 4) phierror*=scalefactors_phi[2];
441  // z error. All clusters have a z-size of 1.
442  const float zerror = length * invsqrt12;
443 
444  double cluslocaly = NAN;
445  double cluslocalz = NAN;
446 
447  if (_make_e_weights[layer])
448  {
449  cluslocaly = ylocalsum / (double) clus_adc;
450  cluslocalz = zlocalsum / (double) clus_adc;
451  }
452  else
453  {
454  cluslocaly = ylocalsum / nhits;
455  cluslocalz = zlocalsum / nhits;
456  }
457 
458  // Fill the cluster fields
459  clus->setAdc(clus_adc);
460 
461  if(Verbosity() > 10) clus->identify();
462 
463  clus->setLocalX(cluslocaly);
464  clus->setLocalY(cluslocalz);
467  clus->setSubSurfKey(0);
468  clus->setActsLocalError(0,0, square(phierror));
469  clus->setActsLocalError(0,1, 0.);
470  clus->setActsLocalError(1,0, 0.);
471  clus->setActsLocalError(1,1, square(zerror));
472  m_clusterlist->addCluster(clus.release());
473 
474  } // end loop over cluster ID's
475  } // end loop over hitsets
476 
477 
478  if(Verbosity() > 1)
479  {
480  // check that the associations were written correctly
481  cout << "After InttClusterizer, cluster-hit associations are:" << endl;
483  }
484 
485  return;
486 }
487 
489 {
490  if (Verbosity() >= 1)
491  {
492  TrkrClusterContainer *clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
493  if (!clusterlist) return;
494 
495  cout << "================= After InttClusterizer::process_event() ====================" << endl;
496 
497  cout << " There are " << clusterlist->size() << " clusters recorded: " << endl;
498 
499  clusterlist->identify();
500 
501  cout << "===========================================================================" << endl;
502  }
503 
504  return;
505 }