31 template<
class T>
inline constexpr
T square(
const T&
x ) {
return x*
x; }
37 if( phi >=
M_PI )
return phi - 2*
M_PI;
38 else if( phi < -M_PI )
return phi + 2*
M_PI;
54 static constexpr
float m_phimin = 0;
55 static constexpr
float m_phimax = 2.*
M_PI;
59 static constexpr
float m_rmin = 20;
60 static constexpr
float m_rmax = 78;
63 static constexpr
float m_zmin = -105.5;
64 static constexpr
float m_zmax = 105.5;
109 <<
"TpcSpaceChargeReconstruction::InitRun\n"
115 <<
" m_max_dz: " <<
m_max_dz <<
"\n"
143 std::unique_ptr<TFile> outputfile( TFile::Open(
m_outputfile.c_str(),
"RECREATE" ) );
150 <<
"TpcSpaceChargeReconstruction::End -"
157 <<
"TpcSpaceChargeReconstruction::End -"
180 m_track_map = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxTrackMap");
183 m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
208 if(
pt < 0.5 )
return false;
211 if( get_clusters<TrkrDefs::mvtxId>(track) < 2 )
return false;
212 if( get_clusters<TrkrDefs::inttId>(track) < 2 )
return false;
213 if(
m_use_micromegas && get_clusters<TrkrDefs::micromegasId>(track) < 2 )
return false;
230 const auto& cluster_key = *key_iter;
234 std::cout <<
PHWHERE <<
" unable to find cluster for key " << cluster_key << std::endl;
245 const auto cluster_r =
get_r( cluster->getX(), cluster->getY() );
246 const auto cluster_phi = std::atan2( cluster->getY(), cluster->getX() );
247 const auto cluster_z = cluster->getZ();
250 const auto cluster_rphi_error = cluster->getRPhiError();
251 const auto cluster_z_error = cluster->getZError();
258 if( cluster_rphi_error < 0.015 )
continue;
259 if( cluster_z_error < 0.05 )
continue;
264 for(
auto iter = state_iter; iter != track->
end_states(); ++iter )
266 const auto dr =
std::abs( cluster_r -
get_r( iter->second->get_x(), iter->second->get_y() ) );
267 if( dr_min < 0 || dr < dr_min )
275 const auto state = state_iter->second;
278 const auto track_r =
get_r( state->get_x(), state->get_y() );
279 const auto dr = cluster_r - track_r;
280 const auto track_drdt = (state->get_x()*state->get_px() + state->get_y()*state->get_py())/track_r;
281 const auto track_dxdr = state->get_px()/track_drdt;
282 const auto track_dydr = state->get_py()/track_drdt;
283 const auto track_dzdr = state->get_pz()/track_drdt;
286 const auto track_x = state->get_x() + dr*track_dxdr;
287 const auto track_y = state->get_y() + dr*track_dydr;
288 const auto track_z = state->get_z() + dr*track_dzdr;
289 const auto track_phi = std::atan2( track_y, track_x );
292 const auto cosphi( std::cos( track_phi ) );
293 const auto sinphi( std::sin( track_phi ) );
294 const auto track_pphi = -state->get_px()*sinphi + state->get_py()*cosphi;
295 const auto track_pr = state->get_px()*cosphi + state->get_py()*sinphi;
296 const auto track_pz = state->get_pz();
297 const auto talpha = -track_pphi/track_pr;
298 const auto tbeta = -track_pz/track_pr;
301 if( std::isnan(talpha) )
303 std::cout <<
"TpcSpaceChargeReconstruction::process_track - talpha is nan" << std::endl;
307 if( std::isnan(tbeta) )
309 std::cout <<
"TpcSpaceChargeReconstruction::process_track - tbeta is nan" << std::endl;
318 const auto track_rphi_error = state->get_rphi_error();
319 const auto track_z_error = state->get_z_error();
322 const auto drp = cluster_r*
delta_phi( cluster_phi - track_phi );
323 const auto dz = cluster_z - track_z;
326 if( std::isnan(drp) )
328 std::cout <<
"TpcSpaceChargeReconstruction::process_track - drp is nan" << std::endl;
334 std::cout <<
"TpcSpaceChargeReconstruction::process_track - dz is nan" << std::endl;
343 const auto erp =
square(track_rphi_error) +
square(cluster_rphi_error);
344 const auto ez =
square(track_z_error) +
square(cluster_z_error);
348 if( std::isnan( erp ) )
350 std::cout <<
"TpcSpaceChargeReconstruction::process_track - erp is nan" << std::endl;
354 if( std::isnan( ez ) )
356 std::cout <<
"TpcSpaceChargeReconstruction::process_track - ez is nan" << std::endl;
364 std::cout <<
"TpcSpaceChargeReconstruction::process_track - invalid cell index" << std::endl;
407 float phi = std::atan2( cluster->
getY(), cluster->
getX() );
408 while( phi < m_phimin ) phi += 2.*
M_PI;
409 while( phi >= m_phimax ) phi -= 2.*
M_PI;
410 int iphi = phibins*(phi-m_phimin)/(m_phimax-m_phimin);
414 if( r < m_rmin || r >= m_rmax )
return -1;
415 int ir = rbins*(r-m_rmin)/(m_rmax-m_rmin);
418 const float z = cluster->
getZ();
419 if( z < m_zmin || z >= m_zmax )
return -1;
420 int iz = zbins*(z-m_zmin)/(m_zmax-m_zmin);