g4simulation/g4eval/TrkrClusterIsMatcher.cc

0001 #include "TrkrClusterIsMatcher.h"
0002 #include "g4evalfn.h"
0003
0004 #include <g4detectors/PHG4CylinderGeomContainer.h>
0005 #include <g4detectors/PHG4TpcCylinderGeom.h>
0006 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0007
0008 #include <intt/CylinderGeomIntt.h>
0009 #include <mvtx/CylinderGeom_Mvtx.h>
0010
0011 #include <trackbase/ActsGeometry.h>
0012 #include <trackbase/TrkrCluster.h>
0013 #include <trackbase/TrkrClusterContainer.h>
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016 #include <phool/PHCompositeNode.h>
0017 #include <phool/PHNode.h>  // for PHNode
0018 #include <phool/PHNodeIterator.h>
0019 #include <phool/PHObject.h>  // for PHObject
0020 #include <phool/getClass.h>
0021 #include <phool/phool.h>  // for PHWHERE
0022
0023 int TrkrClusterIsMatcher::init(
0024     PHCompositeNode* topNode,
0025     const std::string& name_phg4_clusters,
0026     const std::string& name_reco_clusters)
0027 {
0028   // ------ MVTX data ------
0029   PHG4CylinderGeomContainer* geom_container_mvtx = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
0030   if (!geom_container_mvtx)
0031   {
0032     std::cout << PHWHERE << " Could not locate CYLINDERGEOM_MVTX " << std::endl;
0033     return Fun4AllReturnCodes::ABORTRUN;
0034   }
0035   for (int i = 0; i < 3; ++i)
0036   {
0037     auto layergeom = dynamic_cast<CylinderGeom_Mvtx*>(geom_container_mvtx->GetLayerGeom(i));
0038     pitch_phi[i] = layergeom->get_pixel_x();
0039     if (i == 0)
0040     {
0041       pitch_z_MVTX = layergeom->get_pixel_z();
0042     }
0043   }
0044
0045   // ------ INTT data ------
0046   PHG4CylinderGeomContainer* geom_container_intt = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
0047   if (!geom_container_intt)
0048   {
0049     std::cout << PHWHERE << " Could not locate CYLINDERGEOM_INTT " << std::endl;
0050     return Fun4AllReturnCodes::ABORTRUN;
0051   }
0052   // get phi and Z steps for intt
0053   for (int i = 3; i < 7; ++i)
0054   {
0055     CylinderGeomIntt* geom =
0056         dynamic_cast<CylinderGeomIntt*>(geom_container_intt->GetLayerGeom(i));
0057     pitch_phi[i] = geom->get_strip_y_spacing();
0058   }
0059
0060   // ------ TPC data ------
0061   auto geom_tpc =
0062       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0063   if (!geom_tpc)
0064   {
0065     std::cout << PHWHERE << " Could not locate CYLINDERCELLGEOM_SVTX node " << std::endl;
0066     return Fun4AllReturnCodes::ABORTRUN;
0067   }
0068   for (int i = 7; i < 55; ++i)
0069   {
0070     PHG4TpcCylinderGeom* layergeom = geom_tpc->GetLayerCellGeom(i);
0071     if (i == 7)
0072     {
0073       step_t_TPC = layergeom->get_zstep();
0074     }
0075     pitch_phi[i] = layergeom->get_phistep() * layergeom->get_radius();
0076   }
0077
0078   m_TruthClusters =
0079       findNode::getClass<TrkrClusterContainer>(topNode, name_phg4_clusters.c_str());
0080   if (!m_TruthClusters)
0081   {
0082     std::cout << PHWHERE << " Could not locate " << name_phg4_clusters << " node" << std::endl;
0083     return Fun4AllReturnCodes::ABORTRUN;
0084   }
0085
0086   m_RecoClusters =
0087       findNode::getClass<TrkrClusterContainer>(topNode, name_reco_clusters.c_str());
0088   if (!m_TruthClusters)
0089   {
0090     std::cout << PHWHERE << " Could not locate " << name_reco_clusters << " node" << std::endl;
0091     return Fun4AllReturnCodes::ABORTRUN;
0092   }
0093
0094   m_ActsGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0095   if (!m_ActsGeometry)
0096   {
0097     std::cout << PHWHERE << " Could not locate ActsGeometry node" << std::endl;
0098     return Fun4AllReturnCodes::ABORTRUN;
0099   }
0100
0101   set_tol_z_MVTX(tol_z_MVTX);
0102   set_tol_phi_MVTX(tol_phi_MVTX);
0103
0104   set_tol_phi_INTT(tol_phi_INTT);
0105
0106   set_tol_t_TPC(tol_t_TPC);
0107   set_tol_phi_TPC(tol_phi_TPC);
0108
0109   return Fun4AllReturnCodes::EVENT_OK;
0110 }
0111
0112 void TrkrClusterIsMatcher::set_tol_phi_MVTX(float _val)
0113 {
0114   tol_phi_MVTX = _val;
0115   for (int i = 0; i < 3; ++i)
0116   {
0117     tol_pitch_phi[i] = tol_phi_MVTX * pitch_phi[i];
0118   }
0119 }
0120
0121 void TrkrClusterIsMatcher::set_tol_z_MVTX(float _val)
0122 {
0123   tol_z_MVTX = _val;
0124   tol_pitch_z_MVTX = tol_z_MVTX * pitch_z_MVTX;
0125 }
0126
0127 void TrkrClusterIsMatcher::set_tol_phi_INTT(float _val)
0128 {
0129   tol_phi_INTT = _val;
0130   for (int i = 3; i < 7; ++i)
0131   {
0132     tol_pitch_phi[i] = tol_phi_INTT * pitch_phi[i];
0133   }
0134 }
0135
0136 void TrkrClusterIsMatcher::set_tol_phi_TPC(float _val)
0137 {
0138   tol_phi_TPC = _val;
0139   for (int i = 7; i < 55; ++i)
0140   {
0141     tol_pitch_phi[i] = tol_phi_TPC * pitch_phi[i];
0142   }
0143 }
0144
0145 void TrkrClusterIsMatcher::set_tol_t_TPC(float _val)
0146 {
0147   tol_t_TPC = _val;
0148   tol_step_t_TPC = tol_t_TPC * step_t_TPC;
0149 }
0150
0151 bool TrkrClusterIsMatcher::operator()(TrkrDefs::cluskey key_T, TrkrDefs::cluskey key_R)
0152 {
0153   // note: can use returned values, or just pull from these values
0154   layer = TrkrDefs::getLayer(key_T);
0155   if (layer > 55)
0156   {
0157     std::cout << " Error! Trying to compar cluster in layer > 55, "
0158               << "which is not programmed yet!" << std::endl;
0159     return false;
0160   }
0161
0162   clus_T = m_TruthClusters->findCluster(key_T);
0163   clus_R = m_RecoClusters->findCluster(key_R);
0164
0165   det_0123 = g4evalfn::trklayer_det(layer);
0166   dphi = g4evalfn::abs_dphi(clus_T->getPosition(0), clus_R->getPosition(0));
0167   switch (det_0123)
0168   {
0169   case g4evalfn::DET::MVTX:
0170   {
0171     if (single_pixel_phi_MVTX)
0172     {
0173       if (dphi > tol_pitch_phi[layer])
0174       {
0175         return false;
0176       }
0177     }
0178     else
0179     {
0180       if (dphi > tol_pitch_phi[layer] * std::max(clus_T->getPhiSize(), clus_R->getPhiSize()))
0181       {
0182         return false;
0183       }
0184     }
0185     const float delta_z = fabs(clus_T->getPosition(1) - clus_R->getPosition(1));
0186     if (single_pixel_z_MVTX)
0187     {
0188       if (delta_z > tol_pitch_z_MVTX)
0189       {
0190         return false;
0191       }
0192     }
0193     else
0194     {
0195       if (delta_z > tol_pitch_z_MVTX * std::max(clus_T->getZSize(), clus_R->getZSize()))
0196       {
0197         return false;
0198       }
0199     }
0200     return true;
0201   }
0202   break;
0203
0204   case g4evalfn::DET::INTT:
0205   {
0206     if (single_pixel_phi_INTT)
0207     {
0208       if (dphi > tol_pitch_phi[layer])
0209       {
0210         return false;
0211       }
0212     }
0213     else
0214     {
0215       if (dphi > tol_pitch_phi[layer] * std::max(clus_T->getPhiSize(), clus_R->getPhiSize()))
0216       {
0217         return false;
0218       }
0219     }
0220     return true;
0221   }
0222   break;
0223
0224   case g4evalfn::DET::TPC:
0225   {
0226     if (single_bin_phi_TPC)
0227     {
0228       if (dphi > tol_pitch_phi[layer])
0229       {
0230         return false;
0231       }
0232     }
0233     else
0234     {
0235       if (dphi > tol_pitch_phi[layer] * std::max(clus_T->getPhiSize(), clus_R->getPhiSize()))
0236       {
0237         return false;
0238       }
0239     }
0240     const float delta_t = fabs(clus_T->getPosition(1) - clus_R->getPosition(1));
0241     if (single_bin_t_TPC)
0242     {
0243       if (delta_t > tol_step_t_TPC)
0244       {
0245         return false;
0246       }
0247     }
0248     else
0249     {
0250       if (delta_t > tol_step_t_TPC * std::max(clus_T->getZSize(), clus_R->getZSize()))
0251       {
0252         return false;
0253       }
0254     }
0255     return true;
0256   }
0257   break;
0258
0259   case g4evalfn::DET::TPOT:
0260   {                // TPOT
0261     return false;  // no info for matching TPOT at this time
0262   }
0263   break;
0264   }
0265   return false;  // code shouldn't arrive here; just for completeness for compiler
0266 }