packages/trackreco/PHCosmicSeedCombiner.cc

0001
0002 #include "PHCosmicSeedCombiner.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/getClass.h>
0011 #include <phool/phool.h>
0012
0013 #include <trackbase_historic/TrackSeed.h>
0014 #include <trackbase_historic/TrackSeedContainer.h>
0015 #include <trackbase_historic/TrackSeedContainer_v1.h>
0016 #include <trackbase_historic/TrackSeedHelper.h>
0017
0018 #include <trackbase/ActsGeometry.h>
0019 #include <trackbase/TrkrCluster.h>
0020 #include <trackbase/TrkrClusterContainer.h>
0021
0022 #include <cmath>
0023
0024 #include <phool/PHCompositeNode.h>
0025
0026 //____________________________________________________________________________..
0027 PHCosmicSeedCombiner::PHCosmicSeedCombiner(const std::string& name)
0028   : SubsysReco(name)
0029 {
0030 }
0031
0032 //____________________________________________________________________________..
0033 int PHCosmicSeedCombiner::Init(PHCompositeNode*)
0034 {
0035   return Fun4AllReturnCodes::EVENT_OK;
0036 }
0037
0038 //____________________________________________________________________________..
0039 int PHCosmicSeedCombiner::InitRun(PHCompositeNode* topNode)
0040 {
0041   return getNodes(topNode);
0042 }
0043
0044 //____________________________________________________________________________..
0045 int PHCosmicSeedCombiner::process_event(PHCompositeNode*)
0046 {
0047   for (auto trackiter = m_seedMap->begin(); trackiter != m_seedMap->end();
0048        ++trackiter)
0049   {
0050     TrackSeed* track1 = *trackiter;
0051     if (!track1)
0052     {
0053       continue;
0054     }
0055
0056     unsigned int tpcid1 = track1->get_tpc_seed_index();
0057     unsigned int siid1 = track1->get_silicon_seed_index();
0058
0059     if (Verbosity() > 1)
0060     {
0061       std::cout << "tpc id" << tpcid1 << std::endl;
0062     }
0063
0064     auto tpcseed1 = m_tpcSeeds->get(tpcid1);
0065     auto silseed1 = m_siliconSeeds->get(siid1);
0066
0067     const float phi1 = tpcseed1->get_phi();
0068     const float eta1 = tpcseed1->get_eta();
0069
0070     for (auto trackiter2 = trackiter; trackiter2 != m_seedMap->end(); ++trackiter2)
0071     {
0072       if (trackiter == trackiter2) continue;
0073       TrackSeed* track2 = *trackiter2;
0074       if (!track2)
0075       {
0076         continue;
0077       }
0078
0079       unsigned int tpcid2 = track2->get_tpc_seed_index();
0080       unsigned int siid2 = track2->get_silicon_seed_index();
0081
0082       if (Verbosity() > 1)
0083       {
0084         std::cout << "tpc 2 " << tpcid2 << std::endl;
0085       }
0086
0087       auto tpcseed2 = m_tpcSeeds->get(tpcid2);
0088       auto silseed2 = m_siliconSeeds->get(siid2);
0089       const float phi2 = tpcseed2->get_phi();
0090       const float eta2 = tpcseed2->get_eta();
0091
0092       //! phis are the same, so we subtract to check
0093       float dphi = phi1 - phi2;
0094       if (dphi > M_PI)
0095       {
0096         dphi -= 2. * M_PI;
0097       }
0098       else if (dphi < -1 * M_PI)
0099       {
0100         dphi += 2. * M_PI;
0101       }
0102       //! If they are back to back, dphi=pi
0103       dphi = std::abs(dphi) - M_PI;
0104
0105       //! etas are opposite each other, so we add them
0106       const float deta = eta1 + eta2;
0107       if (Verbosity() > 3)
0108       {
0109         std::cout << "phi 1 and phi2  " << phi1 << " , " << phi2 << std::endl;
0110         std::cout << "eta 1 and eta2 " << eta1 << " , " << eta2 << std::endl;
0111         std::cout << "dphi and deta " << dphi << " , " << deta << std::endl;
0112       }
0113       if (std::abs(dphi) < m_dphiCut && std::abs(deta) < m_detaCut)
0114       {
0115         //! add the clusters to the tpc seed and delete seed 2 since it is
0116         //! from the same track
0117         addKeys(tpcseed1, tpcseed2);
0118         if (silseed1)
0119         {
0120           if (silseed2)
0121           {
0122             addKeys(silseed1, silseed2);
0123           }
0124         }
0125         else
0126         {
0127           if (silseed2)
0128           {
0129             track1->set_silicon_seed_index(siid2);
0130           }
0131         }
0132
0133         m_seedMap->erase(m_seedMap->index(trackiter2));
0134       }
0135       if (Verbosity() > 3)
0136       {
0137         track1->identify();
0138         tpcseed1->identify();
0139         if (silseed1)
0140         {
0141           silseed1->identify();
0142         }
0143       }
0144     }
0145   }
0146
0147   return Fun4AllReturnCodes::EVENT_OK;
0148 }
0149 void PHCosmicSeedCombiner::addKeys(TrackSeed* seedToAddTo, TrackSeed* seedToAdd)
0150 {
0151   for (auto citer = seedToAdd->begin_cluster_keys();
0152        citer != seedToAdd->end_cluster_keys();
0153        ++citer)
0154   {
0155     seedToAddTo->insert_cluster_key(*citer);
0156   }
0157 }
0158 //____________________________________________________________________________..
0159 int PHCosmicSeedCombiner::End(PHCompositeNode*)
0160 {
0161   return Fun4AllReturnCodes::EVENT_OK;
0162 }
0163 int PHCosmicSeedCombiner::getNodes(PHCompositeNode* topNode)
0164 {
0165   m_tpcSeeds = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0166   if (!m_tpcSeeds)
0167   {
0168     std::cout << PHWHERE << "TpcTrackSeedContainer not on node tree. Bailing"
0169               << std::endl;
0170     return Fun4AllReturnCodes::ABORTEVENT;
0171   }
0172
0173   m_siliconSeeds = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0174   if (!m_siliconSeeds)
0175   {
0176     std::cout << PHWHERE << "SiliconTrackSeedContainer not on node tree. Bailing"
0177               << std::endl;
0178     return Fun4AllReturnCodes::ABORTEVENT;
0179   }
0180
0181   m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0182   if (!m_clusterContainer)
0183   {
0184     std::cout << PHWHERE
0185               << "No trkr cluster container, exiting." << std::endl;
0186     return Fun4AllReturnCodes::ABORTEVENT;
0187   }
0188
0189   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0190   if (!m_tGeometry)
0191   {
0192     std::cout << "ActsGeometry not on node tree. Exiting."
0193               << std::endl;
0194
0195     return Fun4AllReturnCodes::ABORTEVENT;
0196   }
0197
0198   m_seedMap = findNode::getClass<TrackSeedContainer>(topNode, "SvtxTrackSeedContainer");
0199   if (!m_seedMap)
0200   {
0201     std::cout << "No Svtx seed map on node tree. Exiting."
0202               << std::endl;
0203     return Fun4AllReturnCodes::ABORTEVENT;
0204   }
0205
0206   return Fun4AllReturnCodes::EVENT_OK;
0207 }