ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
MicromegasClusterizer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file MicromegasClusterizer.cc
1 
7 #include "MicromegasDefs.h"
9 
11 #include <g4detectors/PHG4CylinderGeom.h> // for PHG4CylinderGeom
12 
13 #include <trackbase/TrkrClusterContainerv3.h> // for TrkrCluster
15 #include <trackbase/TrkrDefs.h>
16 #include <trackbase/TrkrHitSet.h>
17 #include <trackbase/TrkrHit.h>
20 
21 #include <Acts/Utilities/Units.hpp>
23 
25 #include <fun4all/SubsysReco.h> // for SubsysReco
26 
27 #include <phool/getClass.h>
28 #include <phool/PHCompositeNode.h>
29 #include <phool/PHIODataNode.h> // for PHIODataNode
30 #include <phool/PHNode.h> // for PHNode
31 #include <phool/PHNodeIterator.h> // for PHNodeIterator
32 #include <phool/PHObject.h> // for PHObject
33 
34 #include <Eigen/Dense>
35 
36 #include <TVector3.h>
37 
38 #include <cassert>
39 #include <cmath>
40 #include <cstdint> // for uint16_t
41 #include <iterator> // for distance
42 #include <map> // for _Rb_tree_const_it...
43 #include <utility> // for pair, make_pair
44 #include <vector>
45 
46 
47 namespace
48 {
50  template<class T>
51  inline constexpr T square( const T& x ) { return x*x; }
52 
53  // streamers
54  [[maybe_unused]] inline std::ostream& operator << (std::ostream& out, const Acts::Vector3D& vector )
55  {
56  out << "( " << vector[0] << "," << vector[1] << "," << vector[2] << ")";
57  return out;
58  }
59 
60  // streamers
61  [[maybe_unused]] inline std::ostream& operator << (std::ostream& out, const Acts::Vector2D& vector )
62  {
63  out << "( " << vector[0] << "," << vector[1] << ")";
64  return out;
65  }
66 
67  // streamers
68  [[maybe_unused]] inline std::ostream& operator << (std::ostream& out, const TVector3& vector )
69  {
70  out << "( " << vector.x() << "," << vector.y() << "," << vector.z() << ")";
71  return out;
72  }
73 
74 }
75 
76 //_______________________________________________________________________________
77 MicromegasClusterizer::MicromegasClusterizer(const std::string &name, const std::string& detector)
78  : SubsysReco(name)
79  , m_detector( detector )
80 {}
81 
82 //_______________________________________________________________________________
84 {
85  PHNodeIterator iter(topNode);
86 
87  // Looking for the DST node
88  auto dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
89  assert( dstNode );
90 
91  // Create the Cluster node if missing
92  auto trkrClusterContainer = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
93  if (!trkrClusterContainer)
94  {
95  PHNodeIterator dstiter(dstNode);
96  auto trkrNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
97  if(!trkrNode)
98  {
99  trkrNode = new PHCompositeNode("TRKR");
100  dstNode->addNode(trkrNode);
101  }
102 
103  trkrClusterContainer = new TrkrClusterContainerv3;
104  auto TrkrClusterContainerNode = new PHIODataNode<PHObject>(trkrClusterContainer, "TRKR_CLUSTER", "PHObject");
105  trkrNode->addNode(TrkrClusterContainerNode);
106  }
107 
108  // create cluster to hit association node, if missing
109  auto trkrClusterHitAssoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,"TRKR_CLUSTERHITASSOC");
110  if(!trkrClusterHitAssoc)
111  {
112  PHNodeIterator dstiter(dstNode);
113  auto trkrNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
114  if(!trkrNode)
115  {
116  trkrNode = new PHCompositeNode("TRKR");
117  dstNode->addNode(trkrNode);
118  }
119 
120  trkrClusterHitAssoc = new TrkrClusterHitAssocv3;
121  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(trkrClusterHitAssoc, "TRKR_CLUSTERHITASSOC", "PHObject");
122  trkrNode->addNode(newNode);
123  }
125 }
126 
127 //_______________________________________________________________________________
129 {
130 
131  // geometry
132  PHG4CylinderGeomContainer* geonode = nullptr;
133  for( const std::string& geonodename: {"CYLINDERGEOM_" + m_detector + "_FULL", "CYLINDERGEOM_" + m_detector } )
134  { if(( geonode = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename.c_str()) )) break; }
135  assert(geonode);
136 
137  // hitset container
138  auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
139  assert( trkrhitsetcontainer );
140 
141  // cluster container
142  auto trkrClusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
143  assert( trkrClusterContainer );
144 
145  // cluster-hit association
146  auto trkrClusterHitAssoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
147  assert( trkrClusterHitAssoc );
148 
149  // geometry
150  auto acts_geometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
151  assert( acts_geometry );
152 
153  // surface map
154  auto acts_surface_map = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
155  assert( acts_surface_map );
156 
157  // loop over micromegas hitsets
158  const auto hitset_range = trkrhitsetcontainer->getHitSets(TrkrDefs::TrkrId::micromegasId);
159  for( auto hitset_it = hitset_range.first; hitset_it != hitset_range.second; ++hitset_it )
160  {
161  // get hitset, key and layer
162  TrkrHitSet* hitset = hitset_it->second;
163  const TrkrDefs::hitsetkey hitsetkey = hitset_it->first;
164  const auto layer = TrkrDefs::getLayer(hitsetkey);
165  const auto tileid = MicromegasDefs::getTileId(hitsetkey);
166 
167  // get micromegas geometry object
168  const auto layergeom = dynamic_cast<CylinderGeomMicromegas*>(geonode->GetLayerGeom(layer));
169  assert(layergeom);
170 
171  // get micromegas acts surface
172  const auto acts_surface_iter = acts_surface_map->mmSurfaceMap.find( hitsetkey );
173  if( acts_surface_iter == acts_surface_map->mmSurfaceMap.end() )
174  {
175  std::cout
176  << "MicromegasClusterizer::process_event -"
177  << " could not find surface for layer " << (int) layer << " tile: " << (int) tileid
178  << " skipping hitset"
179  << std::endl;
180  continue;
181  }
182 
183  // surface, surface center and normal director
184  const auto acts_surface( acts_surface_iter->second );
185  Acts::Vector3D normal = acts_surface->normal(acts_geometry->geoContext);
186 
187  if( Verbosity() )
188  {
189  const auto geo_normal = layergeom->get_world_from_local_vect( tileid, {0, 1, 0} );
190  std::cout << "MicromegasClusterizer::process_event -"
191  << " layer: " << (int) layer
192  << " tile: " << (int) tileid
193  << " normal (acts): " << normal
194  << " geo: " << geo_normal
195  << std::endl;
196  }
197 
198  /*
199  * get segmentation type, layer thickness, strip length and pitch.
200  * They are used to calculate cluster errors
201  */
202  const auto segmentation_type = layergeom->get_segmentation_type();
203  const double thickness = layergeom->get_thickness();
204  const double pitch = layergeom->get_pitch();
205  const double strip_length = layergeom->get_strip_length( tileid );
206 
207  // keep a list of ranges corresponding to each cluster
208  using range_list_t = std::vector<TrkrHitSet::ConstRange>;
209  range_list_t ranges;
210 
211  // loop over hits
212  const auto hit_range = hitset->getHits();
213 
214  // keep track of first iterator of runing cluster
215  auto begin = hit_range.first;
216 
217  // keep track of previous strip
218  uint16_t previous_strip = 0;
219  bool first = true;
220 
221  for( auto hit_it = hit_range.first; hit_it != hit_range.second; ++hit_it )
222  {
223 
224  // get hit key
225  const auto hitkey = hit_it->first;
226 
227  // get strip number
228  const auto strip = MicromegasDefs::getStrip( hitkey );
229 
230  if( first )
231  {
232 
233  previous_strip = strip;
234  first = false;
235  continue;
236 
237  } else if( strip - previous_strip > 1 ) {
238 
239  // store current cluster range
240  ranges.push_back( std::make_pair( begin, hit_it ) );
241 
242  // reinitialize begin of next cluster range
243  begin = hit_it;
244 
245  }
246 
247  // update previous strip
248  previous_strip = strip;
249 
250  }
251 
252  // store last cluster
253  if( begin != hit_range.second ) ranges.push_back( std::make_pair( begin, hit_range.second ) );
254 
255  // initialize cluster count
256  int cluster_count = 0;
257 
258  // loop over found hit ranges and create clusters
259  for( const auto& range : ranges )
260  {
261 
262  // create cluster key and corresponding cluster
263  const auto cluster_key = MicromegasDefs::genClusterKey( hitsetkey, cluster_count++ );
264  auto cluster = std::make_unique<TrkrClusterv3>();
265  cluster->setClusKey(cluster_key);
266 
267  TVector3 local_coordinates;
268  double weight_sum = 0;
269 
270  // needed for proper error calculation
271  // it is either the sum over z, or phi, depending on segmentation
272  double coord_sum = 0;
273  double coordsquare_sum = 0;
274 
275  // also store adc value
276  unsigned int adc_sum = 0;
277 
278  // loop over constituting hits
279  for( auto hit_it = range.first; hit_it != range.second; ++hit_it )
280  {
281  // get hit key
282  const auto hitkey = hit_it->first;
283  const auto hit = hit_it->second;
284 
285  // associate cluster key to hit key
286  trkrClusterHitAssoc->addAssoc(cluster_key, hitkey );
287 
288  // get strip number
289  const auto strip = MicromegasDefs::getStrip( hitkey );
290 
291  // get adc, remove pedestal
292  /* pedestal should be the same as the one used in PHG4MicromegasDigitizer */
293  static constexpr double pedestal = 74.6;
294  const double weight = double(hit->getAdc()) - pedestal;
295 
296  // increment cluster adc
297  adc_sum += hit->getAdc();
298 
299  // get strip local coordinate and update relevant sums
300  const auto strip_local_coordinate = layergeom->get_local_coordinates( tileid, strip );
301  local_coordinates += strip_local_coordinate*weight;
302  switch( segmentation_type )
303  {
305  {
306 
307  coord_sum += strip_local_coordinate.x()*weight;
308  coordsquare_sum += square(strip_local_coordinate.x())*weight;
309  break;
310  }
311 
313  {
314  coord_sum += strip_local_coordinate.z()*weight;
315  coordsquare_sum += square(strip_local_coordinate.z())*weight;
316  break;
317  }
318  }
319 
320  weight_sum += weight;
321 
322  }
323 
324  local_coordinates *= (1./weight_sum);
325  const auto world_coordinates = layergeom->get_world_from_local_coords( tileid, local_coordinates);
326  cluster->setAdc( adc_sum );
327 
328  // dimension and error in r, rphi and z coordinates
329  static const float invsqrt12 = 1./std::sqrt(12);
330  static constexpr float error_scale_phi = 1.6;
331  static constexpr float error_scale_z = 0.8;
332 
333  using matrix_t = Eigen::Matrix<float, 3, 3>;
334  matrix_t error = matrix_t::Zero();
335 
336  auto coord_cov = coordsquare_sum/weight_sum - square( coord_sum/weight_sum );
337  auto coord_error_sq = coord_cov/weight_sum;
338  switch( segmentation_type )
339  {
341  if( coord_error_sq == 0 ) coord_error_sq = square(pitch)/12;
342  else coord_error_sq *= square(error_scale_phi);
343  error(0,0) = square(thickness*invsqrt12);
344  error(1,1) = coord_error_sq;
345  error(2,2) = square(strip_length*invsqrt12);
346  break;
347 
349  if( coord_error_sq == 0 ) coord_error_sq = square(pitch)/12;
350  else coord_error_sq *= square(error_scale_z);
351  error(0,0) = square(thickness*invsqrt12);
352  error(1,1) = square(strip_length*invsqrt12);
353  error(2,2) = coord_error_sq;
354  break;
355  }
356 
358  cluster->setLocalX(-1*local_coordinates[0]);
359  cluster->setLocalY(local_coordinates[2]);
360 
361  cluster->setActsLocalError(0,0, error(1,1));
362  cluster->setActsLocalError(0,1, error(1,2));
363  cluster->setActsLocalError(1,0, error(2,1));
364  cluster->setActsLocalError(1,1,error(2,2));
365 
366  // add to container
367  trkrClusterContainer->addCluster( cluster.release() );
368 
369  }
370 
371  }
372  // done
374 }