51 template<
class T>
class range_adaptor
54 range_adaptor(
const T& range ):m_range(range){}
55 inline const typename T::first_type& begin() {
return m_range.first;}
56 inline const typename T::second_type& end() {
return m_range.second;}
62 template<
class T>
inline constexpr
T square(
const T&
x ) {
return x*
x; }
68 template<
class T>
T get_pt(
T px,
T py ) {
return std::sqrt(
square(px) +
square(py) ); }
71 template<
class T>
T get_p(
T px,
T py,
T pz ) {
return std::sqrt(
square(px) +
square(py) +
square(pz) ); }
74 template<
class T>
T get_eta(
T p,
T pz ) {
return std::log( (p+pz)/(p-pz) )/2; }
77 struct interpolation_data_t
79 using list = std::vector<interpolation_data_t>;
84 double px()
const {
return momentum.x(); }
85 double py()
const {
return momentum.y(); }
86 double pz()
const {
return momentum.z(); }
94 template<
double (
interpolation_data_t::*accessor)() const>
95 double average(
const interpolation_data_t::list&
hits )
103 for(
const auto& hit:hits )
106 const double x = (hit.*accessor)();
107 if(std::isnan(x))
continue;
109 const double w = hit.weight;
110 if( w <= 0 )
continue;
117 if( !
valid )
return NAN;
148 trackStruct.
mask = get_mask( track );
157 trackStruct.
x = track->
get_x();
158 trackStruct.
y = track->
get_y();
159 trackStruct.
z = track->
get_z();
160 trackStruct.
r =
get_r( trackStruct.
x, trackStruct.
y );
161 trackStruct.
phi = std::atan2( trackStruct.
y, trackStruct.
x );
166 trackStruct.
pt = get_pt( trackStruct.
px, trackStruct.
py );
167 trackStruct.
p = get_p( trackStruct.
px, trackStruct.
py, trackStruct.
pz );
168 trackStruct.
eta = get_eta( trackStruct.
p, trackStruct.
pz );
176 if( !cluster_hit_map )
return;
177 const auto range = cluster_hit_map->
getHits(clus_key);
180 cluster.
size = std::distance( range.first, range.second );
196 for(
const auto& [first, hit_key]:range_adaptor(range))
222 cluster.
z_size = zbins.size();
234 if(!(cluster_hit_map && hitsetcontainer))
return;
241 const auto hitset = hitsetcontainer->
findHitSet( hitset_key );
242 if( !hitset )
return;
244 const auto range = cluster_hit_map->
getHits(clus_key);
248 for(
const auto& pair:range_adaptor(range))
250 const auto hit = hitset->getHit( pair.second );
253 const auto energy = hit->getEnergy();
278 double line_circle_intersection(
const TVector3& p0,
const TVector3& p1,
double radius )
280 const double A =
square(p1.x() - p0.x()) +
square(p1.y() - p0.y());
281 const double B = 2*p0.x()*(p1.x()-p0.x()) + 2*p0.y()*(p1.y()-p0.y());
284 if( delta < 0 )
return -1;
287 const double tup = (-B + std::sqrt(delta))/(2*A);
288 if( tup >= 0 && tup < 1 )
return tup;
291 const double tdn = (-B-sqrt(delta))/(2*A);
292 if( tdn >= 0 && tdn < 1 )
return tdn;
301 out <<
"ClusterStruct" << std::endl;
302 out <<
" cluster: (" << cluster.
x <<
"," << cluster.
y <<
"," << cluster.
z <<
")" << std::endl;
303 out <<
" track: (" << cluster.
trk_x <<
"," << cluster.
trk_y <<
"," << cluster.
trk_z <<
")" << std::endl;
304 out <<
" truth: (" << cluster.
truth_x <<
"," << cluster.
truth_y <<
"," << cluster.
truth_z <<
")" << std::endl;
310 out <<
"(" << position.x() <<
", " << position.y() <<
", " << position.z() <<
")";
329 std::cout <<
"TrackEvaluation::Init - DST Node missing" << std::endl;
339 std::cout <<
"TrackEvaluation::Init - EVAL node missing - creating" << std::endl;
341 dstNode->addNode(evalNode);
345 evalNode->addNode(newNode);
382 m_surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,
"ActsSurfaceMaps");
386 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
390 m_track_map = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxTrackMap");
393 m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"CORRECTED_TRKR_CLUSTER");
395 { m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER"); }
398 m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
401 m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
404 m_container = findNode::getClass<TrackEvaluationContainerv1>(topNode,
"TrackEvaluationContainer");
407 m_hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
410 m_g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_TPC");
411 m_g4hits_intt = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_INTT");
412 m_g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_MVTX");
416 m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
419 m_tpc_geom_container = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode,
"CYLINDERCELLGEOM_SVTX");
420 assert( m_tpc_geom_container );
450 const int nclusters = std::distance( clusters.first, clusters.second );
460 event.nclusters[
layer] += nclusters;
517 const auto track = trackpair.second;
518 auto track_struct = create_track( track );
522 track_struct.contributors = contributors;
526 track_struct.embed =
get_embed(particle);
530 auto state_iter = track->begin_states();
533 for(
auto key_iter = track->begin_cluster_keys(); key_iter != track->end_cluster_keys(); ++key_iter )
536 const auto& cluster_key = *key_iter;
540 std::cout <<
"TrackEvaluation::evaluate_tracks - unable to find cluster for key " << cluster_key << std::endl;
562 const auto radius( cluster_struct.r );
564 for(
auto iter = state_iter; iter != track->end_states(); ++iter )
566 const auto dr =
std::abs( radius -
get_r( iter->second->get_x(), iter->second->get_y() ) );
567 if( dr_min < 0 || dr < dr_min )
584 track_struct.clusters.push_back( cluster_struct );
598 auto map_iter =
m_g4hit_map.lower_bound( cluster_key );
599 if( map_iter !=
m_g4hit_map.end() && cluster_key == map_iter->first )
600 {
return map_iter->second; }
613 for(
auto truth_iter = g4hit_map.begin(); truth_iter != g4hit_map.end(); ++truth_iter )
616 const auto g4hit_key = truth_iter->second.second;
640 if( g4hit ) out.insert( g4hit );
641 else std::cout <<
"TrackEvaluation::find_g4hits - g4hit not found " << g4hit_key << std::endl;
647 return m_g4hit_map.insert( map_iter, std::make_pair( cluster_key, std::move( out ) ) )->second;
657 using IdMap = std::map<int,int>;
658 IdMap contributor_map;
663 const auto& cluster_key = *key_iter;
666 const int trkid = hit->get_trkid();
667 auto iter = contributor_map.lower_bound( trkid );
668 if( iter == contributor_map.end() || iter->first != trkid )
670 contributor_map.insert(iter, std::make_pair(trkid,1));
671 }
else ++iter->second;
675 if( contributor_map.empty() )
return {0,0};
676 else return *std::max_element(
677 contributor_map.cbegin(), contributor_map.cend(),
678 [](
const IdMap::value_type& first,
const IdMap::value_type&
second )
679 {
return first.second <
second.second; } );
695 cluster_struct.
x = global.x();
696 cluster_struct.
y = global.y();
697 cluster_struct.
z = global.z();
698 cluster_struct.
r =
get_r( cluster_struct.
x, cluster_struct.
y );
699 cluster_struct.
phi = std::atan2( cluster_struct.
y, cluster_struct.
x );
703 return cluster_struct;
711 const auto dr = cluster.
r - trk_r;
713 const auto trk_dxdr = state->
get_px()/trk_drdt;
714 const auto trk_dydr = state->
get_py()/trk_drdt;
715 const auto trk_dzdr = state->
get_pz()/trk_drdt;
733 const auto cosphi( std::cos( cluster.
trk_phi ) );
734 const auto sinphi( std::sin( cluster.
trk_phi ) );
735 const auto trk_pphi = -state->
get_px()*sinphi + state->
get_py()*cosphi;
736 const auto trk_pr = state->
get_px()*cosphi + state->
get_py()*sinphi;
737 const auto trk_pz = state->
get_pz();
738 cluster.
trk_alpha = std::atan2( trk_pphi, trk_pr );
739 cluster.
trk_beta = std::atan2( trk_pz, trk_pr );
755 const TVector3 cluster_world( cluster.
x, cluster.
y, cluster.
z );
756 const TVector3 cluster_local = layergeom->get_local_from_world_coords( tileid, cluster_world );
759 TVector3 track_world( state->
get_x(), state->
get_y(), state->
get_z() );
760 TVector3 track_local = layergeom->get_local_from_world_coords( tileid, track_world );
764 const TVector3 direction_local = layergeom->get_local_from_world_vect( tileid, direction_world );
767 const auto delta_y = cluster_local.y() - track_local.y();
768 track_local += TVector3(
769 delta_y*direction_local.x()/direction_local.y(),
771 delta_y*direction_local.z()/direction_local.y() );
774 track_world = layergeom->get_world_from_local_coords( tileid, track_local );
777 cluster.
trk_x = track_world.x();
778 cluster.
trk_y = track_world.y();
779 cluster.
trk_z = track_world.z();
792 const auto cosphi( std::cos( cluster.
trk_phi ) );
793 const auto sinphi( std::sin( cluster.
trk_phi ) );
794 const auto trk_pphi = -state->
get_px()*sinphi + state->
get_py()*cosphi;
795 const auto trk_pr = state->
get_px()*cosphi + state->
get_py()*sinphi;
796 const auto trk_pz = state->
get_pz();
797 cluster.
trk_alpha = std::atan2( trk_pphi, trk_pr );
798 cluster.
trk_beta = std::atan2( trk_pz, trk_pr );
812 const bool is_tpc(
layer >= 7 &&
layer < 55 );
818 interpolation_data_t::list
hits;
819 for(
const auto& g4hit:g4hits )
821 interpolation_data_t::list tmp_hits;
822 const auto weight = g4hit->get_edep();
823 for(
int i = 0; i < 2; ++i )
825 const TVector3 g4hit_world(g4hit->get_x(i), g4hit->get_y(i), g4hit->get_z(i));
826 const TVector3 momentum_world(g4hit->get_px(i), g4hit->get_py(i), g4hit->get_pz(i));
827 tmp_hits.push_back( {.position = g4hit_world, .momentum = momentum_world, .weight =
weight } );
834 auto r0 =
get_r(tmp_hits[0].
x(),tmp_hits[0].
y());
835 auto r1 =
get_r(tmp_hits[1].
x(),tmp_hits[1].
y());
843 if( r1 <= rin || r0 >= rout )
continue;
846 const auto dr_old =
r1-r0;
851 const auto t = line_circle_intersection( tmp_hits[0].position, tmp_hits[1].position, rin );
854 tmp_hits[0].position = tmp_hits[0].position*(1.-
t) + tmp_hits[1].position*
t;
855 tmp_hits[0].momentum = tmp_hits[0].momentum*(1.-
t) + tmp_hits[1].
momentum*
t;
861 const auto t = line_circle_intersection( tmp_hits[0].position, tmp_hits[1].position, rout );
864 tmp_hits[1].position = tmp_hits[0].position*(1.-
t) + tmp_hits[1].position*
t;
865 tmp_hits[1].momentum = tmp_hits[0].momentum*(1.-
t) + tmp_hits[1].
momentum*
t;
870 const auto dr_new =
r1-r0;
871 tmp_hits[0].weight *= dr_new/dr_old;
872 tmp_hits[1].weight *= dr_new/dr_old;
876 hits.push_back(std::move(tmp_hits[0]));
877 hits.push_back(std::move(tmp_hits[1]));
881 cluster.
truth_x = average<&interpolation_data_t::x>(
hits );
882 cluster.
truth_y = average<&interpolation_data_t::y>(
hits );
883 cluster.
truth_z = average<&interpolation_data_t::z>(
hits );
886 cluster.
truth_px = average<&interpolation_data_t::px>(
hits );
887 cluster.
truth_py = average<&interpolation_data_t::py>(
hits );
888 cluster.
truth_pz = average<&interpolation_data_t::pz>(
hits );
897 const auto cosphi( std::cos( cluster.
truth_phi ) );
898 const auto sinphi( std::sin( cluster.
truth_phi ) );
902 cluster.
truth_alpha = std::atan2( truth_pphi, truth_pr );
918 const TVector3 cluster_world( cluster.
x, cluster.
y, cluster.
z );
919 const TVector3 cluster_local = layergeom->get_local_from_world_coords( tileid, cluster_world );
922 interpolation_data_t::list
hits;
923 for(
const auto& g4hit:g4hits )
925 const auto weight = g4hit->get_edep();
929 layergeom->convert_to_planar( tileid, &g4hit_copy );
931 for(
int i = 0; i < 2; ++i )
935 TVector3 g4hit_world(g4hit_copy.
get_x(i), g4hit_copy.
get_y(i), g4hit_copy.
get_z(i));
936 TVector3 g4hit_local = layergeom->get_local_from_world_coords( tileid, g4hit_world );
939 TVector3 momentum_world(g4hit_copy.
get_px(i), g4hit_copy.
get_py(i), g4hit_copy.
get_pz(i));
940 TVector3 momentum_local = layergeom->get_local_from_world_vect( tileid, momentum_world );
942 hits.push_back( {.position = g4hit_local, .momentum = momentum_local, .weight =
weight } );
947 const TVector3 interpolation_local(
948 average<&interpolation_data_t::x>(hits),
949 average<&interpolation_data_t::y>(hits),
950 average<&interpolation_data_t::z>(hits) );
952 const TVector3 interpolation_world = layergeom->get_world_from_local_coords( tileid, interpolation_local );
955 const TVector3 momentum_local(
956 average<&interpolation_data_t::px>(hits),
957 average<&interpolation_data_t::py>(hits),
958 average<&interpolation_data_t::pz>(hits));
960 const TVector3 momentum_world = layergeom->get_world_from_local_vect( tileid, momentum_local );
962 cluster.
truth_x = interpolation_world.x();
963 cluster.
truth_y = interpolation_world.y();
964 cluster.
truth_z = interpolation_world.z();
969 cluster.
truth_px = momentum_world.x();
970 cluster.
truth_py = momentum_world.y();
971 cluster.
truth_pz = momentum_world.z();
977 const auto cosphi( std::cos( cluster.
truth_phi ) );
978 const auto sinphi( std::sin( cluster.
truth_phi ) );
982 cluster.
truth_alpha = std::atan2( truth_pphi, truth_pr );