g4simulation/g4eval/TrkrClusterIsMatcher.cc

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