34 #include <Eigen/Dense>
51 inline constexpr
T square(
const T&
x ) {
return x*
x; }
56 out <<
"( " << vector[0] <<
"," << vector[1] <<
"," << vector[2] <<
")";
63 out <<
"( " << vector[0] <<
"," << vector[1] <<
")";
68 [[maybe_unused]]
inline std::ostream&
operator << (std::ostream& out,
const TVector3& vector )
70 out <<
"( " << vector.x() <<
"," << vector.y() <<
"," << vector.z() <<
")";
79 , m_detector( detector )
92 auto trkrClusterContainer = findNode::getClass<TrkrClusterContainer>(dstNode,
"TRKR_CLUSTER");
93 if (!trkrClusterContainer)
100 dstNode->addNode(trkrNode);
104 auto TrkrClusterContainerNode =
new PHIODataNode<PHObject>(trkrClusterContainer,
"TRKR_CLUSTER",
"PHObject");
105 trkrNode->addNode(TrkrClusterContainerNode);
109 auto trkrClusterHitAssoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
110 if(!trkrClusterHitAssoc)
117 dstNode->addNode(trkrNode);
122 trkrNode->addNode(newNode);
133 for(
const std::string& geonodename: {
"CYLINDERGEOM_" +
m_detector +
"_FULL",
"CYLINDERGEOM_" +
m_detector } )
134 {
if(( geonode = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename.c_str()) ))
break; }
138 auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
139 assert( trkrhitsetcontainer );
142 auto trkrClusterContainer = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
143 assert( trkrClusterContainer );
146 auto trkrClusterHitAssoc = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
147 assert( trkrClusterHitAssoc );
150 auto acts_geometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
151 assert( acts_geometry );
154 auto acts_surface_map = findNode::getClass<ActsSurfaceMaps>(topNode,
"ActsSurfaceMaps");
155 assert( acts_surface_map );
159 for(
auto hitset_it = hitset_range.first; hitset_it != hitset_range.second; ++hitset_it )
172 const auto acts_surface_iter = acts_surface_map->mmSurfaceMap.find( hitsetkey );
173 if( acts_surface_iter == acts_surface_map->mmSurfaceMap.end() )
176 <<
"MicromegasClusterizer::process_event -"
177 <<
" could not find surface for layer " << (
int)
layer <<
" tile: " << (
int) tileid
178 <<
" skipping hitset"
184 const auto acts_surface( acts_surface_iter->second );
189 const auto geo_normal = layergeom->get_world_from_local_vect( tileid, {0, 1, 0} );
190 std::cout <<
"MicromegasClusterizer::process_event -"
192 <<
" tile: " << (
int) tileid
193 <<
" normal (acts): " << normal
194 <<
" geo: " << geo_normal
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 );
208 using range_list_t = std::vector<TrkrHitSet::ConstRange>;
212 const auto hit_range = hitset->
getHits();
215 auto begin = hit_range.first;
218 uint16_t previous_strip = 0;
221 for(
auto hit_it = hit_range.first; hit_it != hit_range.second; ++hit_it )
225 const auto hitkey = hit_it->first;
233 previous_strip = strip;
237 }
else if( strip - previous_strip > 1 ) {
240 ranges.push_back( std::make_pair( begin, hit_it ) );
248 previous_strip = strip;
253 if( begin != hit_range.second ) ranges.push_back( std::make_pair( begin, hit_range.second ) );
256 int cluster_count = 0;
259 for(
const auto& range : ranges )
264 auto cluster = std::make_unique<TrkrClusterv3>();
265 cluster->setClusKey(cluster_key);
267 TVector3 local_coordinates;
268 double weight_sum = 0;
272 double coord_sum = 0;
273 double coordsquare_sum = 0;
276 unsigned int adc_sum = 0;
279 for(
auto hit_it = range.first; hit_it != range.second; ++hit_it )
282 const auto hitkey = hit_it->first;
283 const auto hit = hit_it->second;
286 trkrClusterHitAssoc->addAssoc(cluster_key,
hitkey );
293 static constexpr
double pedestal = 74.6;
294 const double weight = double(hit->getAdc()) - pedestal;
297 adc_sum += hit->getAdc();
300 const auto strip_local_coordinate = layergeom->get_local_coordinates( tileid, strip );
301 local_coordinates += strip_local_coordinate*
weight;
302 switch( segmentation_type )
307 coord_sum += strip_local_coordinate.x()*
weight;
308 coordsquare_sum +=
square(strip_local_coordinate.x())*weight;
314 coord_sum += strip_local_coordinate.z()*
weight;
315 coordsquare_sum +=
square(strip_local_coordinate.z())*weight;
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 );
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;
333 using matrix_t = Eigen::Matrix<float, 3, 3>;
334 matrix_t
error = matrix_t::Zero();
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 )
341 if( coord_error_sq == 0 ) coord_error_sq =
square(pitch)/12;
342 else coord_error_sq *=
square(error_scale_phi);
344 error(1,1) = coord_error_sq;
349 if( coord_error_sq == 0 ) coord_error_sq =
square(pitch)/12;
350 else coord_error_sq *=
square(error_scale_z);
353 error(2,2) = coord_error_sq;
358 cluster->setLocalX(-1*local_coordinates[0]);
359 cluster->setLocalY(local_coordinates[2]);
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));
367 trkrClusterContainer->addCluster( cluster.release() );