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(
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; }
152 trackStruct.
mask = get_mask( track );
161 trackStruct.
x = track->
get_x();
162 trackStruct.
y = track->
get_y();
163 trackStruct.
z = track->
get_z();
164 trackStruct.
r =
get_r( trackStruct.
x, trackStruct.
y );
165 trackStruct.
phi = std::atan2( trackStruct.
y, trackStruct.
x );
170 trackStruct.
pt = get_pt( trackStruct.
px, trackStruct.
py );
171 trackStruct.
p = get_p( trackStruct.
px, trackStruct.
py, trackStruct.
pz );
172 trackStruct.
eta = get_eta( trackStruct.
p, trackStruct.
pz );
182 cluster_struct.
x = cluster->
getX();
183 cluster_struct.
y = cluster->
getY();
184 cluster_struct.
z = cluster->
getZ();
185 cluster_struct.
r =
get_r( cluster_struct.
x, cluster_struct.
y );
186 cluster_struct.
phi = std::atan2( cluster_struct.
y, cluster_struct.
x );
189 std::cout <<
" (x|y|z|r|l) "
190 << cluster_struct.
x <<
" | "
191 << cluster_struct.
y <<
" | "
192 << cluster_struct.
z <<
" | "
193 << cluster_struct.
r <<
" | "
194 << cluster_struct.
layer <<
" | "
196 std::cout <<
" (xl|yl) "
200 return cluster_struct;
208 const auto dr = cluster.
r - trk_r;
210 const auto trk_dxdr = state->
get_px()/trk_drdt;
211 const auto trk_dydr = state->
get_py()/trk_drdt;
212 const auto trk_dzdr = state->
get_pz()/trk_drdt;
230 const auto cosphi( std::cos( cluster.
trk_phi ) );
231 const auto sinphi( std::sin( cluster.
trk_phi ) );
232 const auto trk_pphi = -state->
get_px()*sinphi + state->
get_py()*cosphi;
233 const auto trk_pr = state->
get_px()*cosphi + state->
get_py()*sinphi;
234 const auto trk_pz = state->
get_pz();
235 cluster.
trk_alpha = std::atan2( trk_pphi, trk_pr );
236 cluster.
trk_beta = std::atan2( trk_pz, trk_pr );
245 if( !cluster_hit_map )
return;
246 const auto range = cluster_hit_map->
getHits(clus_key);
249 cluster.
size = std::distance( range.first, range.second );
265 for(
const auto& [first, hit_key]:range_adaptor(range))
291 cluster.
z_size = zbins.size();
303 if(!(cluster_hit_map && hitsetcontainer))
return;
310 const auto hitset = hitsetcontainer->
findHitSet( hitset_key );
311 if( !hitset )
return;
313 const auto range = cluster_hit_map->
getHits(clus_key);
317 for(
const auto& pair:range_adaptor(range))
319 const auto hit = hitset->getHit( pair.second );
322 const auto energy = hit->getEnergy();
335 const auto rextrap = cluster.
r;
337 std::cout <<
"inter" <<
"\n";
339 cluster.
truth_x = interpolate<&PHG4Hit::get_x>(
hits, rextrap );
340 cluster.
truth_y = interpolate<&PHG4Hit::get_y>(
hits, rextrap );
341 cluster.
truth_z = interpolate<&PHG4Hit::get_z>(
hits, rextrap );
346 cluster.
truth_px = interpolate<&PHG4Hit::get_px>(
hits, rextrap );
347 cluster.
truth_py = interpolate<&PHG4Hit::get_py>(
hits, rextrap );
348 cluster.
truth_pz = interpolate<&PHG4Hit::get_pz>(
hits, rextrap );
350 std::cout <<
"inter2" <<
"\n";
356 const auto cosphi( std::cos( cluster.
truth_phi ) );
357 const auto sinphi( std::sin( cluster.
truth_phi ) );
361 cluster.
truth_alpha = std::atan2( truth_pphi, truth_pr );
366 double truth_alpha = 0;
367 double truth_beta = 0;
369 for(
const auto& hit:hits )
371 const auto px = hit->get_x(1) - hit->get_x(0);
372 const auto py = hit->get_y(1) - hit->get_y(0);
373 const auto pz = hit->get_z(1) - hit->get_z(0);
374 const auto pphi = -px*sinphi + py*cosphi;
375 const auto pr = px*cosphi + py*sinphi;
377 const auto w = hit->get_edep();
378 if(
w < 0 )
continue;
381 truth_alpha +=
w*std::atan2( pphi,
pr );
382 truth_beta +=
w*std::atan2( pz,
pr );
384 truth_alpha /= sum_w;
414 , _filename(filename)
417 , sabotage(inSabotage)
418 , apply_compression(compress)
427 _dst_data =
new TNtuple(
"dst_data",
"dst data",
"event:seed:"
428 "c3x:c3y:c3z:c3p:t3x:t3y:t3z:"
429 "c2x:c2y:c2r:c2l:t2x:t2y:t2r:t2l:"
432 "pt:eta:phi:charge");
439 std::cout <<
"DSTEmulator::Init - DST Node missing" << std::endl;
449 std::cout <<
"DSTEmulator::Init - EVAL node missing - creating" << std::endl;
451 dstNode->addNode(evalNode);
455 evalNode->addNode(newNode);
487 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
488 "ActsTrackingGeometry");
492 <<
"ActsTrackingGeometry not found on node tree. Exiting"
497 m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
502 <<
"ActsSurfaceMaps not found on node tree. Exiting"
545 m_track_map = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxTrackMapTPCOnly");
548 std::cout <<
" DSTEmulator: Using TPC Only Track Map node " << std::endl;
552 std::cout <<
" DSTEmulator: TPC Only Track Map node not found, using default" << std::endl;
553 m_track_map = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxTrackMap");
557 m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"CORRECTED_TRKR_CLUSTER");
561 std::cout <<
" DSTEmulator: Using CORRECTED_TRKR_CLUSTER node " << std::endl;
565 std::cout <<
" DSTEmulator: CORRECTED_TRKR_CLUSTER node not found, using TRKR_CLUSTER" << std::endl;
566 m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
570 std::cout <<
" DSTEmulator: CORRECTED_TRKR_CLUSTER node not found at all, using TRKR_CLUSTER" << std::endl;
571 m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
576 m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
579 m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
582 m_container = findNode::getClass<TrackEvaluationContainerv1>(topNode,
"TrackEvaluationContainer");
585 m_hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
588 m_g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_TPC");
589 m_g4hits_intt = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_INTT");
590 m_g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_MVTX");
594 m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
618 const auto track = trackpair.second;
619 auto track_struct = create_track( track );
623 track_struct.contributors = contributors;
627 track_struct.embed =
get_embed(particle);
631 auto state_iter = track->begin_states();
634 for(
auto key_iter = track->begin_cluster_keys(); key_iter != track->end_cluster_keys(); ++key_iter )
637 const auto& cluster_key = *key_iter;
643 std::cout <<
"DSTEmulator::evaluate_tracks - unable to find cluster for key " << cluster_key << std::endl;
656 float radius =
get_r( globalpos_d[0], globalpos_d[1] );
657 float clu_phi = std::atan2( globalpos_d[0], globalpos_d[1] );
658 std::cout <<
"radius " << radius << std::endl;
660 for(
auto iter = state_iter; iter != track->end_states(); ++iter )
662 const auto dr =
std::abs( radius -
get_r( iter->second->get_x(), iter->second->get_y() ) );
663 if( dr_min < 0 || dr < dr_min )
673 std::cout <<
"NEW (xg|yg) "
674 << globalpos_d[0] <<
" | "
677 std::cout <<
"NEW (xl|yl) "
678 << cluster->getLocalX() <<
" | "
679 << cluster->getLocalY()
687 const auto trk_r =
get_r( state_iter->second->get_x(), state_iter->second->get_y() );
688 std::cout <<
" trk_r " << trk_r << std::endl;
689 const auto dr =
get_r( globalpos_d[0], globalpos_d[1] ) - trk_r;
690 std::cout <<
" dr " << dr << std::endl;
691 const auto trk_drdt = (state_iter->second->get_x()*state_iter->second->get_px() + state_iter->second->get_y()*state_iter->second->get_py())/trk_r;
692 std::cout <<
" trk_drdt " << trk_drdt << std::endl;
693 const auto trk_dxdr = state_iter->second->get_px()/trk_drdt;
694 std::cout <<
" trk_dxdr " << trk_dxdr << std::endl;
695 const auto trk_dydr = state_iter->second->get_py()/trk_drdt;
696 std::cout <<
" trk_dydr " << trk_dydr << std::endl;
697 const auto trk_dzdr = state_iter->second->get_pz()/trk_drdt;
698 std::cout <<
" trk_dzdr " << trk_dzdr << std::endl;
702 float trk_x = state_iter->second->get_x() + dr*trk_dxdr;
703 float trk_y = state_iter->second->get_y() + dr*trk_dydr;
704 float trk_z = state_iter->second->get_z() + dr*trk_dzdr;
705 std::cout <<
"trk_x " << state_iter->second->get_x() <<
"trk_y" << state_iter->second->get_y() <<
"trk_z " << state_iter->second->get_z() << std::endl;
738 std::map<unsigned int, std::vector<Surface>>::iterator mapIter;
743 <<
"Error: hitsetkey not found in clusterSurfaceMap, layer = " <<
746 << hitsetkey << std::endl;
749 std::cout <<
" g0: " << global[0] <<
" g1: " << global[1] <<
" g2:" << global[2] << std::endl;
756 std::vector<Surface> surf_vec = mapIter->second;
759 double world_phi = atan2(world[1], world[0]);
760 double world_z = world[2];
762 double fraction = (world_phi +
M_PI) / (2.0 *
M_PI);
763 double rounded_nsurf =
round( (
double) (surf_vec.size()/2) * fraction - 0.5);
764 unsigned int nsurf = (
unsigned int) rounded_nsurf;
766 nsurf += surf_vec.size()/2;
775 double TrkRadius = sqrt(trk_x * trk_x + trk_y * trk_y);
776 double rTrkPhi = TrkRadius * atan2(trk_y, trk_x);
777 double surfRadius = sqrt(center(0)*center(0) + center(1)*center(1));
778 double surfPhiCenter = atan2(center[1], center[0]);
779 double surfRphiCenter = surfPhiCenter * surfRadius;
780 double surfZCenter = center[2];
786 delta_r = TrkRadius - surfRadius;
787 trklocX = rTrkPhi - surfRphiCenter;
788 trklocY = trk_z - surfZCenter;
800 float delta_x = trklocX - cluster->getLocalX();
801 float delta_y = trklocY - cluster->getLocalY();
808 comp_dx =
rnd.Uniform(-1, 1) * delta_x;
809 comp_dy =
rnd.Uniform(-1, 1) * delta_y;
811 comp_dx =
rnd.Uniform(-10, 10);
812 comp_dy =
rnd.Uniform(-10, 10);
824 (
float)globalpos_d[0],
825 (float)globalpos_d[1],
826 (
float)globalpos_d[2],
831 cluster->getLocalX(),
832 cluster->getLocalY(),
833 get_r((
float)globalpos_d[0],(
float)globalpos_d[1]),
847 (
float) track_struct.charge
850 std::cout <<
"filled" <<
"\n";
853 cluster->setLocalX(trklocX - comp_dx);
854 cluster->setLocalY(trklocY - comp_dy);
878 std::cout <<
"DSTEmulator::evaluate_tracks - tracks: " <<
m_container->
tracks().size() << std::endl;
907 auto map_iter =
m_g4hit_map.lower_bound( cluster_key );
908 if( map_iter !=
m_g4hit_map.end() && cluster_key == map_iter->first )
909 {
return map_iter->second; }
922 for(
auto truth_iter = g4hit_map.begin(); truth_iter != g4hit_map.end(); ++truth_iter )
925 const auto g4hit_key = truth_iter->second.second;
949 if( g4hit ) out.insert( g4hit );
950 else std::cout <<
"DSTEmulator::find_g4hits - g4hit not found " << g4hit_key << std::endl;
956 return m_g4hit_map.insert( map_iter, std::make_pair( cluster_key, std::move( out ) ) )->second;
966 using IdMap = std::map<int,int>;
967 IdMap contributor_map;
972 const auto& cluster_key = *key_iter;
975 const int trkid = hit->get_trkid();
976 auto iter = contributor_map.lower_bound( trkid );
977 if( iter == contributor_map.end() || iter->first != trkid )
979 contributor_map.insert(iter, std::make_pair(trkid,1));
980 }
else ++iter->second;
984 if( contributor_map.empty() )
return {0,0};
985 else return *std::max_element(
986 contributor_map.cbegin(), contributor_map.cend(),
987 [](
const IdMap::value_type& first,
const IdMap::value_type&
second )
988 {
return first.second <
second.second; } );