34 template<
class T>
class range_adaptor
37 range_adaptor(
const T& range ):m_range(range){}
38 inline const typename T::first_type& begin() {
return m_range.first;}
39 inline const typename T::second_type& end() {
return m_range.second;}
46 inline constexpr
T square(
const T&
x ) {
return x*
x; }
53 double line_circle_intersection(
const TVector3&
p,
const TVector3&
d,
double radius )
56 const double B = 2*p.x()*d.x() + 2*p.y()*d.y();
59 if( delta < 0 )
return -1;
62 const double tup = (-B + std::sqrt(delta))/(2*A);
63 if( tup >= 0 )
return tup;
66 const double tdn = (-B-sqrt(delta))/(2*A);
67 if( tdn >= 0 )
return tdn;
74 inline std::ostream&
operator << (std::ostream& out,
const TVector3& vector )
76 out <<
"( " << vector.x() <<
", " << vector.y() <<
", " << vector.z() <<
")";
84 if( phi >=
M_PI )
return phi - 2*
M_PI;
85 else if( phi < -M_PI )
return phi + 2*
M_PI;
90 static constexpr
float m_phimin = 0;
91 static constexpr
float m_phimax = 2.*
M_PI;
95 static constexpr
float m_rmin = 20;
96 static constexpr
float m_rmax = 78;
99 static constexpr
float m_zmin = -105.5;
100 static constexpr
float m_zmax = 105.5;
135 <<
"TpcDirectLaserReconstruction::InitRun\n"
139 <<
" m_max_dz: " <<
m_max_dz <<
"\n"
166 std::unique_ptr<TFile> outputfile( TFile::Open(
m_outputfile.c_str(),
"RECREATE" ) );
176 {
if( o ) o->Write(); }
182 <<
"TpcDirectLaserReconstruction::End -"
215 m_surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,
"ActsSurfaceMaps");
219 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
223 m_track_map = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxTrackMap");
227 m_hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
230 m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
231 assert(m_cluster_map);
238 std::cout <<
"TpcDirectLaserReconstruction::makeHistograms - writing evaluation histograms to: " <<
m_histogramfilename << std::endl;
243 h_dca_layer =
new TH2F(
"dca_layer",
";layer; DCA (cm)", 57, 0, 57, 500, 0, 2 );
244 h_deltarphi_layer =
new TH2F(
"deltarphi_layer",
";layer; r.#Delta#phi_{track-cluster} (cm)", 57, 0, 57, 2000, -2, 2 );
245 h_deltaz_layer =
new TH2F(
"deltaz_layer",
";layer; #Deltaz_{track-cluster} (cm)", 57, 0, 57, 400, -2, 2 );
255 h_entries =
new TH3F(
"entries",
";#phi;r (cm);z (cm)", phibins, m_phimin, m_phimax, rbins, m_rmin, m_rmax, zbins, m_zmin, m_zmax );
287 { std::cout <<
"TpcDirectLaserReconstruction::process_track - position: " <<
origin <<
" direction: " << direction << std::endl; }
307 const TVector3 oc( global.x()-
origin.x(), global.y()-
origin.y(), global.z()-
origin.z() );
308 auto t = direction.Dot( oc )/
square( direction.Mag() );
309 auto om = direction*
t;
310 const auto dca = (oc-
om).Mag();
316 t = line_circle_intersection(
origin, direction,
get_r( global.x(), global.y() ));
317 if( t < 0 )
continue;
323 const auto pathlength = om.Mag();
326 const auto projection =
origin +
om;
330 state.
set_x( projection.x() );
331 state.
set_y( projection.y() );
332 state.
set_z( projection.z() );
334 state.
set_px( direction.x());
335 state.
set_py( direction.y());
336 state.
set_pz( direction.z());
343 const auto cluster_r =
get_r(global.x(), global.y());
344 const auto cluster_phi = std::atan2(global.y(),global.x());
345 const auto cluster_z = global.z();
348 const auto cluster_rphi_error = cluster->getRPhiError();
349 const auto cluster_z_error = cluster->getZError();
360 const auto track_phi = std::atan2( projection.y(), projection.x() );
361 const auto track_z = projection.z();
364 const auto cosphi( std::cos( track_phi ) );
365 const auto sinphi( std::sin( track_phi ) );
366 const auto track_pphi = -state.
get_px()*sinphi + state.
get_py()*cosphi;
367 const auto track_pr = state.
get_px()*cosphi + state.
get_py()*sinphi;
368 const auto track_pz = state.
get_pz();
369 const auto talpha = -track_pphi/track_pr;
370 const auto tbeta = -track_pz/track_pr;
373 if( std::isnan(talpha) )
375 std::cout <<
"TpcDirectLaserReconstruction::process_track - talpha is nan" << std::endl;
379 if( std::isnan(tbeta) )
381 std::cout <<
"TpcDirectLaserReconstruction::process_track - tbeta is nan" << std::endl;
386 const auto drp = cluster_r*
delta_phi( cluster_phi - track_phi );
387 const auto dz = cluster_z - track_z;
390 if( std::isnan(drp) )
392 std::cout <<
"TpcDirectLaserReconstruction::process_track - drp is nan" << std::endl;
398 std::cout <<
"TpcDirectLaserReconstruction::process_track - dz is nan" << std::endl;
409 auto phi = cluster_phi;
410 while( phi < m_phimin ) phi += 2.*
M_PI;
411 while( phi >= m_phimax ) phi -= 2.*
M_PI;
412 h_entries->Fill( phi, cluster_r, cluster_z );
421 const auto erp =
square(cluster_rphi_error);
422 const auto ez =
square(cluster_z_error);
425 if( std::isnan( erp ) )
427 std::cout <<
"TpcDirectLaserReconstruction::process_track - erp is nan" << std::endl;
431 if( std::isnan( ez ) )
433 std::cout <<
"TpcDirectLaserReconstruction::process_track - ez is nan" << std::endl;
443 std::cout <<
"TpcDirectLaserReconstruction::process_track - invalid cell index"
444 <<
" r: " << cluster_r
445 <<
" phi: " << cluster_phi
446 <<
" z: " << cluster_z
493 float phi = std::atan2( global.y(), global.x() );
494 while( phi < m_phimin ) phi += 2.*
M_PI;
495 while( phi >= m_phimax ) phi -= 2.*
M_PI;
496 int iphi = phibins*(phi-m_phimin)/(m_phimax-m_phimin);
499 const float r =
get_r( global.x(), global.y() );
500 if( r < m_rmin || r >= m_rmax )
return -1;
501 int ir = rbins*(r-m_rmin)/(m_rmax-m_rmin);
504 const float z = global.z();
505 if( z < m_zmin || z >= m_zmax )
return -1;
506 int iz = zbins*(z-m_zmin)/(m_zmax-m_zmin);