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