packages/tpc/LaserClusterHelper.cc

0001 #include "LaserClusterHelper.h"
0002
0003 #include <trackbase/LaserCluster.h>
0004 #include <trackbase/TpcDefs.h>
0005 #include <trackbase/TrkrDefs.h>
0006
0007 #include <g4detectors/PHG4TpcGeom.h>
0008 #include <g4detectors/PHG4TpcGeomContainer.h>
0009
0010 #include <phool/PHCompositeNode.h>
0011 #include <phool/getClass.h>
0012
0013 #include <cmath>
0014 #include <iostream>
0015 #include <limits>
0016
0017 namespace
0018 {
0019   Acts::Vector3 invalid(std::numeric_limits<double>::quiet_NaN(),
0020                         std::numeric_limits<double>::quiet_NaN(),
0021                         std::numeric_limits<double>::quiet_NaN());
0022
0023   std::array<double, 3> invalidArr = {std::numeric_limits<double>::quiet_NaN(),
0024                         std::numeric_limits<double>::quiet_NaN(),
0025                         std::numeric_limits<double>::quiet_NaN()};
0026 }
0027
0028 //____________________________________________________________________________
0029 void LaserClusterHelper::loadNodes(PHCompositeNode* topNode)
0030 {
0031     m_tGeometry = findNode::getClass<ActsGeometry>(topNode,"ActsGeometry");
0032     if(!m_tGeometry)
0033     {
0034         std::cout << "LaserClusterHelper::loadNodes - ActsGeometry not found on node tree" << std::endl;
0035     }
0036
0037     m_geom_container = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0038     if(!m_geom_container)
0039     {
0040         std::cout << "LaserClusterHelper::loadNodes - TPCGEOMCONTAINER not found on node tree" << std::endl;
0041     }
0042 }
0043
0044 //____________________________________________________________________________
0045 Acts::Vector3 LaserClusterHelper::getHitPosition(TrkrDefs::hitsetkey hitsetkey, TrkrDefs::hitkey hitkey) const
0046 {
0047     //const Acts::Vector3 invalid(std::numeric_limits<double>::quiet_NaN(),
0048     //                            std::numeric_limits<double>::quiet_NaN(),
0049     //                            std::numeric_limits<double>::quiet_NaN());
0050
0051     if(!m_tGeometry || !m_geom_container)
0052     {
0053         return invalid;
0054     }
0055
0056     const int layer = TrkrDefs::getLayer(hitsetkey);
0057     const int side = TpcDefs::getSide(hitsetkey);
0058
0059     PHG4TpcGeom *layer_geom = m_geom_container->GetLayerCellGeom(layer);
0060     if(!layer_geom)
0061     {
0062         return invalid;
0063     }
0064
0065     const int iphi = TpcDefs::getPad(hitkey);
0066     const int it = TpcDefs::getTBin(hitkey);
0067
0068     const double radius = layer_geom->get_radius();
0069     const double phi = layer_geom->get_phi(iphi, side);
0070
0071     const double env_x = radius * cos(phi);
0072     const double env_y = radius * sin(phi);
0073     double env_z = 0.0;
0074     //hard code at 0 until better z coordinate calibration is determined
0075     if(m_useZ)
0076     {
0077         double vdrift = m_tGeometry->get_drift_velocity();
0078         double tdriftmax = layer_geom->get_max_driftlength() / vdrift;
0079
0080         double zdriftlength = layer_geom->get_zcenter(it) * vdrift;
0081         // convert z drift length to z position in the TPC
0082         env_z = tdriftmax * vdrift - zdriftlength;
0083         if (side == 0)
0084         {
0085             env_z = -env_z;
0086         }
0087     }
0088
0089     Acts::Vector3 env_global(env_x, env_y, env_z);
0090     if(!m_useGlobal)
0091     {
0092         return env_global;
0093     }
0094
0095     return m_tGeometry->transformTpcEnvelopeToWorld(env_global);
0096 }
0097
0098 //____________________________________________________________________________
0099 Acts::Vector3 LaserClusterHelper::getClusterCentroid(LaserCluster* cluster) const
0100 {
0101     //const Acts::Vector3 invalid(std::numeric_limits<double>::quiet_NaN(),
0102     //                            std::numeric_limits<double>::quiet_NaN(),
0103     //                            std::numeric_limits<double>::quiet_NaN());
0104
0105     if(!cluster)
0106     {
0107         return invalid;
0108     }
0109
0110     Acts::Vector3 weightedSum(0.0, 0.0, 0.0);
0111     double adcSum = 0.0;
0112
0113     const unsigned int nhits = cluster->getNhits();
0114     for(unsigned int i=0; i<nhits; ++i)
0115     {
0116         const LaserClusterHitInfo hit= cluster->getHit(i);
0117         const Acts::Vector3 hitCoords = getHitPosition(hit.hitsetkey, hit.hitkey);
0118         if(hitCoords.hasNaN())
0119         {
0120             continue;
0121         }
0122
0123         weightedSum += hit.adc * hitCoords;
0124         adcSum += hit.adc;
0125     }
0126
0127     if(adcSum <= 0.0)
0128     {
0129         return invalid;
0130     }
0131
0132     return weightedSum / adcSum;
0133 }
0134
0135 //____________________________________________________________________________
0136 std::array<double, 3> LaserClusterHelper::getClusterHardwareCentroid(LaserCluster* cluster) const
0137 {
0138     //const Acts::Vector3 invalid(std::numeric_limits<double>::quiet_NaN(),
0139     //                            std::numeric_limits<double>::quiet_NaN(),
0140     //                            std::numeric_limits<double>::quiet_NaN());
0141
0142     if(!cluster)
0143     {
0144         return invalidArr;
0145     }
0146
0147     Acts::Vector3 weightedSum(0.0, 0.0, 0.0);
0148     double adcSum = 0.0;
0149     double layerSum = 0.0;
0150     double iphiSum = 0.0;
0151     double itSum = 0.0;
0152
0153     const unsigned int nhits = cluster->getNhits();
0154     for(unsigned int i=0; i<nhits; ++i)
0155     {
0156         const LaserClusterHitInfo hit= cluster->getHit(i);
0157
0158         const int layer = TrkrDefs::getLayer(hit.hitsetkey);
0159         const int iphi = TpcDefs::getPad(hit.hitkey);
0160         const int it = TpcDefs::getTBin(hit.hitkey);
0161
0162         adcSum += hit.adc;
0163         layerSum += layer * hit.adc;
0164         iphiSum += iphi * hit.adc;
0165         itSum += it * hit.adc;
0166     }
0167
0168     if(adcSum <= 0.0)
0169     {
0170         return invalidArr;
0171     }
0172
0173     return {layerSum / adcSum, iphiSum / adcSum, itSum / adcSum};
0174 }