packages/trackreco/PHTpcDeltaZCorrection.cc

0001 #include "PHTpcDeltaZCorrection.h"
0002
0003 #include "PHTpcDeltaZCorrection.h"
0004
0005 /// Tracking includes
0006
0007 #include <trackbase/TrkrCluster.h>  // for TrkrCluster
0008 #include <trackbase/TrkrClusterContainer.h>
0009 #include <trackbase_historic/TrackSeed.h>
0010 #include <trackbase_historic/TrackSeedContainer.h>
0011 #include <trackbase_historic/TrackSeedHelper.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014
0015 #include <phool/getClass.h>
0016 #include <phool/phool.h>
0017
0018 #include <gsl/gsl_const_mksa.h>  // for the speed of light
0019 #include <cmath>                 // for sqrt, fabs, atan2, cos
0020 #include <iostream>              // for operator<<, basic_ostream
0021 #include <map>                   // for map
0022 #include <set>                   // for _Rb_tree_const_iterator
0023 #include <utility>               // for pair, make_pair
0024
0025 namespace
0026 {
0027   /// square
0028   template <class T>
0029   inline constexpr T square(const T& x)
0030   {
0031     return x * x;
0032   }
0033
0034   /// speed of light, in cm per ns
0035   static constexpr double speed_of_light = GSL_CONST_MKSA_SPEED_OF_LIGHT * 1e-7;
0036
0037   [[maybe_unused]] std::ostream& operator << (std::ostream& out, const Acts::Vector3& v )
0038   {
0039     out << "(" << v.x() << ", " << v.y() << ", " << v.z() << ")";
0040     return out;
0041   }
0042
0043
0044   [[maybe_unused]] std::ostream& operator << (std::ostream& out, const Acts::Vector2& v )
0045   {
0046     out << "(" << v.x() << ", " << v.y() << ")";
0047     return out;
0048   }
0049 }  // namespace
0050
0051 //____________________________________________________________________________..
0052 PHTpcDeltaZCorrection::PHTpcDeltaZCorrection(const std::string& name)
0053   : SubsysReco(name)
0054   , PHParameterInterface(name)
0055 {
0056   InitializeParameters();
0057 }
0058
0059 //____________________________________________________________________________..
0060 int PHTpcDeltaZCorrection::InitRun(PHCompositeNode* /*unused*/)
0061 {
0062   UpdateParametersWithMacro();
0063   return Fun4AllReturnCodes::EVENT_OK;
0064 }
0065
0066 //____________________________________________________________________________..
0067 int PHTpcDeltaZCorrection::process_event(PHCompositeNode* topNode)
0068 {
0069   // load nodes
0070   const auto res = load_nodes(topNode);
0071   if (res != Fun4AllReturnCodes::EVENT_OK)
0072   {
0073     return res;
0074   }
0075
0076   process_tracks();
0077   return Fun4AllReturnCodes::EVENT_OK;
0078 }
0079
0080 //_____________________________________________________________________
0081 int PHTpcDeltaZCorrection::End(PHCompositeNode* /*topNode*/)
0082 {
0083   return Fun4AllReturnCodes::EVENT_OK;
0084 }
0085
0086 //_____________________________________________________________________
0087 void PHTpcDeltaZCorrection::SetDefaultParameters()
0088 {
0089   // Data on gasses @20 C and 760 Torr from the following source:
0090   // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf
0091   // diffusion and drift velocity for 400kV for NeCF4 50/50 from calculations:
0092   // http://skipper.physics.sunysb.edu/~prakhar/tpc/HTML_Gases/split.html
0093   return;
0094 }
0095
0096 //_____________________________________________________________________
0097 int PHTpcDeltaZCorrection::load_nodes(PHCompositeNode* topNode)
0098 {
0099   // acts geometry
0100   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0101   assert(m_tGeometry);
0102
0103   // get necessary nodes
0104   m_track_map = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0105   assert(m_track_map);
0106
0107   m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
0108   if (m_cluster_map)
0109   {
0110     if (Verbosity() > 0)
0111     {
0112       std::cout << " Using CORRECTED_TRKR_CLUSTER node " << std::endl;
0113     }
0114   }
0115   else
0116   {
0117     if (Verbosity() > 0)
0118     {
0119       std::cout << " CORRECTED_TRKR_CLUSTER node not found, using TRKR_CLUSTER" << std::endl;
0120     }
0121     m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, m_clusterContainerName);
0122   }
0123   assert(m_cluster_map);
0124   return Fun4AllReturnCodes::EVENT_OK;
0125 }
0126
0127 //_____________________________________________________________________
0128 void PHTpcDeltaZCorrection::process_tracks()
0129 {
0130   if (!(m_track_map && m_cluster_map))
0131   {
0132     return;
0133   }
0134   for (unsigned int iter = 0; iter != m_track_map->size(); ++iter)
0135   {
0136     TrackSeed* seed = m_track_map->get(iter);
0137     if (!seed)
0138     {
0139       continue;
0140     }
0141     process_track(iter, seed);
0142   }
0143
0144   m_corrected_clusters.clear();
0145 }
0146
0147 //_____________________________________________________________________
0148 void PHTpcDeltaZCorrection::process_track(unsigned int key, TrackSeed* track)
0149 {
0150   // keep track of the global position of previous cluster on track
0151   const auto origin = TrackSeedHelper::get_xyz(track);
0152
0153   // radius
0154   const double radius = fabs(1. / track->get_qOverR());  // cm
0155
0156   // helix center
0157   const double center_x = track->get_X0();
0158   const double center_y = track->get_Y0();
0159
0160   // origin to center 2D vector
0161   const Acts::Vector2 orig_vect = {origin.x() - center_x, origin.y() - center_y};
0162
0163   // print
0164   if (Verbosity())
0165   {
0166     std::cout << "PHTpcDeltaZCorrection -"
0167               << " track: " << key
0168               << " positive: " << track->get_charge()
0169               << " center: " << center_x << ", " << center_y
0170               << " radius: " << radius
0171               << std::endl;
0172   }
0173
0174   /*
0175    * check track seed fit validity
0176    * skip if invalid, otherwise it will result in clusters with NAN coordinates
0177    */
0178   if( std::isnan( center_x ) || std::isnan( center_y ) )
0179   {
0180     std::cout << "PHTpcDeltaZCorrection::process_track - invalid seed parameters. Skipping" << std::endl;
0181     return;
0182   }
0183
0184   // loop over clusters. Assume they are ordered by layer
0185   for (auto key_iter = track->begin_cluster_keys(); key_iter != track->end_cluster_keys(); ++key_iter)
0186   {
0187     // store cluster key
0188     const auto& cluster_key = *key_iter;
0189
0190     // consider TPC clusters only
0191     if (TrkrDefs::getTrkrId(cluster_key) != TrkrDefs::tpcId)
0192     {
0193       continue;
0194     }
0195
0196     // skip if cluster was already corrected
0197     if (!m_corrected_clusters.insert(cluster_key).second)
0198     {
0199       continue;
0200     }
0201
0202     // get cluster
0203     auto cluster = m_cluster_map->findCluster(cluster_key);
0204     if (!cluster)
0205     {
0206       continue;
0207     }
0208
0209     // get cluster global position
0210     const auto global = m_tGeometry->getGlobalPosition(cluster_key, cluster);
0211
0212     // get delta z
0213     const double delta_z = global.z() - origin.z();
0214
0215     // cluster position to center
0216     const Acts::Vector2 cluster_vect = {global.x() - center_x, global.y() - center_y};
0217
0218     // delta phi
0219     const double delta_phi = std::atan2(
0220         cluster_vect.y() * orig_vect.x() - cluster_vect.x() * orig_vect.y(),
0221         cluster_vect.x() * orig_vect.x() + cluster_vect.y() * orig_vect.y());
0222
0223     // helical path length
0224     const double pathlength = std::sqrt(square(delta_z) + square(radius * delta_phi));
0225
0226     // adjust cluster position to account for particles propagation time
0227     /*
0228      * accounting for particles finite velocity results in reducing the electron drift time by pathlenght/c
0229      * this in turn affects the cluster z, so that it is always closer to the readout plane
0230      */
0231     const double t_correction = pathlength / speed_of_light;
0232     cluster->setLocalY(cluster->getLocalY() - t_correction);
0233
0234     if (Verbosity())
0235     {
0236       std::cout << "PHTpcDeltaZCorrection::process_track - cluster: " << cluster_key
0237                 << " path length: " << pathlength << " t correction " << t_correction << std::endl;
0238     }
0239   }
0240 }