File indexing completed on 2025-08-06 08:18:23
0001
0002 #include "AzimuthalSeeder.h"
0003
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <phool/PHCompositeNode.h>
0006 #include <phool/PHDataNode.h>
0007 #include <phool/PHNode.h>
0008 #include <phool/PHNodeIterator.h>
0009 #include <phool/PHObject.h>
0010 #include <phool/PHTimer.h>
0011 #include <phool/getClass.h>
0012 #include <phool/phool.h>
0013
0014 #include <trackbase/TrackFitUtils.h>
0015 #include <trackbase/TrkrCluster.h>
0016 #include <trackbase/TrkrClusterContainer.h>
0017 #include <trackbase_historic/TrackSeed.h>
0018 #include <trackbase_historic/TrackSeedContainer.h>
0019 #include <trackbase_historic/TrackSeedContainer_v1.h>
0020 #include <trackbase_historic/TrackSeed_v2.h>
0021 #include <trackbase_historic/TrackSeedHelper.h>
0022
0023 #include <TFile.h>
0024 #include <TH2.h>
0025
0026 #include <cmath>
0027
0028 AzimuthalSeeder::AzimuthalSeeder(const std::string &name)
0029 : SubsysReco(name)
0030 {
0031 }
0032
0033
0034 int AzimuthalSeeder::Init(PHCompositeNode * )
0035 {
0036 file = std::make_unique<TFile>("ASoutfile.root", "recreate");
0037 h_phi = new TH2F("h_phi", "h_phi", 1000, -3.2, 3.2, 1000, -3.2, 3.2);
0038 h_phi2 = new TH2F("h_phi2", "h_phi2", 1000, -3.2, 3.2, 1000, -3.2, 3.2);
0039 h_phi3 = new TH2F("h_phi3", "h_phi3", 1000, -3.2, 3.2, 1000, -3.2, 3.2);
0040
0041 return Fun4AllReturnCodes::EVENT_OK;
0042 }
0043
0044
0045 int AzimuthalSeeder::InitRun(PHCompositeNode *topNode)
0046 {
0047 int ret = createNodes(topNode);
0048 if (ret != Fun4AllReturnCodes::EVENT_OK)
0049 {
0050 return ret;
0051 }
0052 getNodes(topNode);
0053 return ret;
0054 }
0055
0056
0057 int AzimuthalSeeder::process_event(PHCompositeNode * )
0058 {
0059 AzimuthalSeeder::PositionMap clusterPositions[4];
0060 std::set<TrkrDefs::TrkrId> detectors;
0061 detectors.insert(TrkrDefs::TrkrId::mvtxId);
0062 for (const auto &det : detectors)
0063 {
0064 for (const auto &layer : {0, 1, 2})
0065 {
0066 for (const auto &hitsetkey : m_clusterContainer->getHitSetKeys(det, layer))
0067 {
0068 auto range = m_clusterContainer->getClusters(hitsetkey);
0069 for (auto citer = range.first; citer != range.second; ++citer)
0070 {
0071 const auto ckey = citer->first;
0072 const auto cluster = citer->second;
0073 const auto global = m_tGeometry->getGlobalPosition(ckey, cluster);
0074 clusterPositions[layer].insert(std::make_pair(ckey, global));
0075 }
0076 }
0077 }
0078 }
0079
0080 SeedVector seeds;
0081 for (auto iter = clusterPositions[0].begin(); iter != clusterPositions[0].end(); ++iter)
0082 {
0083 auto key0 = iter->first;
0084 auto pos0 = iter->second;
0085 for (auto &iter2 : clusterPositions[1])
0086 {
0087 auto key1 = iter2.first;
0088 auto pos1 = iter2.second;
0089 if (key0 == key1)
0090 {
0091 continue;
0092 }
0093 float phi0 = atan2(pos0.y(), pos0.x());
0094 float phi1 = atan2(pos1.y(), pos1.x());
0095
0096 h_phi->Fill(phi0, phi0 - phi1);
0097 if (std::fabs(phi0 - phi1) < 0.1)
0098 {
0099 seed s;
0100 s.ckeys.push_back(key0);
0101 s.ckeys.push_back(key1);
0102 s.globpos.push_back(pos0);
0103 s.globpos.push_back(pos1);
0104 seeds.push_back(s);
0105 }
0106 }
0107 }
0108 for (auto &s : seeds)
0109 {
0110 for (const auto &[key2, pos2] : clusterPositions[2])
0111 {
0112 float phi2 = atan2(pos2.y(), pos2.x());
0113 float phi1 = atan2(s.globpos[1].y(), s.globpos[1].x());
0114 float phi0 = atan2(s.globpos[0].y(), s.globpos[0].x());
0115 h_phi2->Fill(phi0, phi0 - phi2);
0116 h_phi3->Fill(phi2, phi1 - phi2);
0117 if (std::fabs(phi0 - phi2) < 0.1 && std::fabs(phi1 - phi2) < 0.1)
0118 {
0119 s.ckeys.push_back(key2);
0120 s.globpos.push_back(pos2);
0121 }
0122 }
0123 }
0124
0125 SeedVector finalseeds;
0126 for (auto &s : seeds)
0127 {
0128 TrackFitUtils::position_vector_t xypoints, rzpoints;
0129 for (auto &pos : s.globpos)
0130 {
0131 xypoints.push_back(std::make_pair(pos.x(), pos.y()));
0132 float r = std::sqrt(pos.x() * pos.x() + pos.y() * pos.y());
0133 if (pos.y() < 0)
0134 {
0135 r = -r;
0136 }
0137 rzpoints.push_back(std::make_pair(pos.z(), r));
0138 }
0139 auto xyLineParams = TrackFitUtils::line_fit(xypoints);
0140 auto rzLineParams = TrackFitUtils::line_fit(rzpoints);
0141 float xySlope = std::get<0>(xyLineParams);
0142 float xyIntercept = std::get<1>(xyLineParams);
0143 float rzSlope = std::get<0>(rzLineParams);
0144 float rzIntercept = std::get<1>(rzLineParams);
0145 bool badseed = false;
0146
0147 for (auto &pos : s.globpos)
0148 {
0149 float r = std::sqrt(pos.x() * pos.x() + pos.y() * pos.y());
0150 if (pos.y() < 0)
0151 {
0152 r = -r;
0153 }
0154 float distancexy = std::abs(xySlope * pos.x() - pos.y() + xyIntercept) / std::sqrt(xySlope * xySlope + 1);
0155 float distancerz = std::abs(rzSlope * pos.z() - r + rzIntercept) / std::sqrt(rzSlope * rzSlope + 1);
0156 if (distancexy > m_outlierLimit || distancerz > m_outlierLimit)
0157 {
0158 badseed = true;
0159 break;
0160 }
0161 }
0162 if (!badseed)
0163 {
0164 finalseeds.push_back(std::move(s));
0165 }
0166 }
0167
0168 for (auto &map : {clusterPositions[1], clusterPositions[2]})
0169 {
0170 for (auto &[key, pos] : map)
0171 {
0172 clusterPositions[0].insert(std::make_pair(key, pos));
0173 }
0174 }
0175 if(Verbosity() > 2)
0176 {
0177 std::cout << "finalseed size " << finalseeds.size() << std::endl;
0178 }
0179 for (auto &s : finalseeds)
0180 {
0181 if (s.ckeys.size() < 3)
0182 {
0183 continue;
0184 }
0185 std::unique_ptr<TrackSeed_v2> si_seed = std::make_unique<TrackSeed_v2>();
0186 for (auto &key : s.ckeys)
0187 {
0188 si_seed->insert_cluster_key(key);
0189 }
0190 float avgphi = 0;
0191 for (auto &pos : s.globpos)
0192 {
0193 avgphi += atan2(pos.y(), pos.x());
0194 }
0195 avgphi /= s.globpos.size();
0196 si_seed->set_phi(avgphi);
0197 TrackSeedHelper::circleFitByTaubin(si_seed.get(),clusterPositions[0], 0, 3);
0198 TrackSeedHelper::lineFit(si_seed.get(),clusterPositions[0], 0, 3);
0199 if (Verbosity() > 3)
0200 {
0201 si_seed->identify();
0202 }
0203 m_seedContainer->insert(si_seed.get());
0204 }
0205 return Fun4AllReturnCodes::EVENT_OK;
0206 }
0207
0208
0209 int AzimuthalSeeder::End(PHCompositeNode * )
0210 {
0211 if (m_outfile)
0212 {
0213 file->cd();
0214 h_phi->Write();
0215 h_phi2->Write();
0216 h_phi3->Write();
0217 file->Close();
0218 }
0219 return Fun4AllReturnCodes::EVENT_OK;
0220 }
0221
0222 int AzimuthalSeeder::createNodes(PHCompositeNode *topNode)
0223 {
0224 PHNodeIterator iter(topNode);
0225
0226 PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0227
0228 if (!dstNode)
0229 {
0230 std::cerr << "DST node is missing, quitting" << std::endl;
0231 throw std::runtime_error("Failed to find DST node in PHSiliconCosmicSeeding::createNodes");
0232 }
0233
0234 PHNodeIterator dstIter(dstNode);
0235 PHCompositeNode *svtxNode = dynamic_cast<PHCompositeNode *>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0236
0237 if (!svtxNode)
0238 {
0239 svtxNode = new PHCompositeNode("SVTX");
0240 dstNode->addNode(svtxNode);
0241 }
0242
0243 m_seedContainer = findNode::getClass<TrackSeedContainer>(topNode, m_trackMapName);
0244 if (!m_seedContainer)
0245 {
0246 m_seedContainer = new TrackSeedContainer_v1;
0247 PHIODataNode<PHObject> *trackNode =
0248 new PHIODataNode<PHObject>(m_seedContainer, m_trackMapName, "PHObject");
0249 svtxNode->addNode(trackNode);
0250 }
0251
0252 return Fun4AllReturnCodes::EVENT_OK;
0253 }
0254
0255 int AzimuthalSeeder::getNodes(PHCompositeNode *topNode)
0256 {
0257 m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0258 if (!m_tGeometry)
0259 {
0260 std::cout << PHWHERE << "No acts reco geometry, bailing."
0261 << std::endl;
0262 return Fun4AllReturnCodes::ABORTEVENT;
0263 }
0264 m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0265 if (!m_clusterContainer)
0266 {
0267 std::cout << PHWHERE << "No cluster container on node tree, bailing."
0268 << std::endl;
0269 return Fun4AllReturnCodes::ABORTEVENT;
0270 }
0271 return Fun4AllReturnCodes::EVENT_OK;
0272 }
