18 #include <calobase/RawTowerGeomContainer.h>
19 #include <calobase/RawTowerContainer.h>
20 #include <calobase/RawTower.h>
21 #include <calobase/RawClusterContainer.h>
22 #include <calobase/RawCluster.h>
23 #include <calobase/RawClusterUtility.h>
24 #include <phgeom/PHGeomUtility.h>
26 #include <Acts/Geometry/GeometryIdentifier.hpp>
33 #include <ActsExamples/Plugins/BField/ScalableBField.hpp>
40 , m_trajectories(nullptr)
54 std::cout <<
"PHActsTrackProjection begin Init" << std::endl;
62 std::cout <<
"PHActsTrackProjection finished Init" << std::endl;
70 std::cout <<
"PHActsTrackProjection : Starting process_event event "
76 std::cout <<
"Processing calo layer "
85 std::cout <<
"PHActsTrackProjection : Finished process_event event "
114 const auto& [trackTips, mj] = traj.trajectory();
115 if(trackTips.size() > 1 and
Verbosity() > 0)
118 <<
"More than 1 track tip per track. Should never happen..."
121 for(
const auto& trackTip : trackTips)
123 const auto params = traj.trackParameters(trackTip);
131 auto trackStateParams = std::move(**result);
151 Acts::Surface::makeShared<Acts::PerigeeSurface>(actsVertex);
159 for(
int i = 0; i < 6; i++)
160 for(
int j = 0; j < 6; j++)
188 auto projectionPhi = atan2(projectionPos(1), projectionPos(0));
189 auto projectionEta = asinh(projectionPos(2) /
190 sqrt(projectionPos(0) * projectionPos(0)
191 + projectionPos(1) * projectionPos(1)));
196 double radius = sqrt(projectionPos(0) * projectionPos(0)
197 + projectionPos(1) * projectionPos(1));
198 std::cout <<
"Track projection phi/eta " << projectionPhi
199 <<
" and " << projectionEta << std::endl
200 <<
" projection position " << projectionPos.transpose()
201 <<
" and radius is " << radius << std::endl;
204 if(fabs(projectionEta) >= 1.1)
210 double energy3x3 = 0.0;
211 double energy5x5 = 0.0;
216 std::cout <<
"3x3/5x5 energy sums " << energy3x3
217 <<
" and " << energy5x5 << std::endl;
224 double minIndex = NAN;
225 double minDphi = NAN;
226 double minDeta = NAN;
230 minIndex, minDphi, minDeta, minE);
232 if(!std::isnan(minIndex))
243 std::cout <<
"Calo cluster has dphi " << minDphi
244 <<
" and deta " << minDeta <<
" and clusID "
245 << minIndex <<
" and energy " << minE
262 for(
const auto& [key, cluster] : clusterMap)
264 const auto clusterEta =
267 const auto dphi =
deltaPhi(phi - cluster->get_phi());
268 const auto deta = eta - clusterEta;
269 const auto r = sqrt(pow(dphi,2) + pow(deta,2));
277 minE = cluster->get_energy();
289 for(
int iphi = phiBin - 2; iphi <= phiBin + 2; ++iphi)
291 for(
int ieta = etaBin - 2; ieta <= etaBin + 2; ++ieta)
310 if(
abs(iphi - phiBin) <= 1 and
abs(ieta - etaBin) <= 1)
311 energy3x3 += tower->get_energy();
325 std::cout <<
"Propagating final track fit with momentum: "
326 << params.momentum() <<
" and position "
329 <<
"track fit phi/eta "
330 << atan2(params.momentum()(1),
331 params.momentum()(0))
333 << atanh(params.momentum()(2)
334 / params.momentum().norm())
337 return std::visit([params, targetSurf,
this]
339 using InputMagneticField =
340 typename std::decay_t<decltype(inputField)>::element_type;
345 MagneticField field(inputField);
358 Acts::LoggerWrapper{*logger});
360 auto result = propagator.propagate(params, *targetSurf,
364 return std::move((*result).endParameters);
366 return result.error();
368 std::move(m_tGeometry->magField));
375 std::string towerGeoNodeName =
"TOWERGEOM_" +
m_caloNames.at(caloLayer);
376 std::string towerNodeName =
"TOWER_CALIB_" +
m_caloNames.at(caloLayer);
377 std::string clusterNodeName =
"CLUSTER_" +
m_caloNames.at(caloLayer);
380 (topNode, towerGeoNodeName.c_str());
383 (topNode, towerNodeName.c_str());
386 (topNode, clusterNodeName.c_str());
391 std::string nodeName =
"CLUSTER_POS_COR_" +
m_caloNames.at(caloLayer);
392 m_clusterContainer = findNode::getClass<RawClusterContainer>
393 (topNode, nodeName.c_str());
399 <<
"Calo geometry and/or cluster container not found on node tree. Bailing."
409 for(
int caloLayer = 0; caloLayer <
m_nCaloLayers; caloLayer++)
418 const auto eta = 2.5;
419 const auto theta = 2. * atan(exp(-
eta));
420 const auto halfZ = caloRadius / tan(
theta)
424 auto transform = Acts::Transform3D::Identity();
426 std::shared_ptr<Acts::CylinderSurface> surf =
427 Acts::Surface::makeShared<Acts::CylinderSurface>(
transform,
439 std::cout <<
"Cylinder " <<
name <<
" has center "
451 m_trajectories = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode,
"ActsTrajectories");
454 std::cout <<
PHWHERE <<
"No Acts trajectories on node tree, bailing."
459 m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode,
"SvtxVertexMap");
462 std::cout <<
PHWHERE <<
"No vertex map on node tree, bailing."
467 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(
468 topNode,
"ActsTrackingGeometry");
471 std::cout <<
"ActsTrackingGeometry not on node tree. Exiting."
477 m_trackMap = findNode::getClass<SvtxTrackMap>(topNode,
"SvtxTrackMap");
480 std::cout <<
PHWHERE <<
"No SvtxTrackMap on node tree. Bailing."
491 return phi - 2. *
M_PI;
492 else if (phi <= -
M_PI)
493 return phi + 2.*
M_PI;