g4simulation/g4eval/SvtxTruthEval.cc

0001 #include "SvtxTruthEval.h"
0002
0003 #include "BaseTruthEval.h"
0004
0005 #include <g4main/PHG4Hit.h>
0006 #include <g4main/PHG4HitContainer.h>
0007 #include <g4main/PHG4TruthInfoContainer.h>
0008
0009 #include <g4main/PHG4Particle.h>
0010
0011 #include <micromegas/MicromegasDefs.h>
0012
0013 #include <trackbase/ActsGeometry.h>
0014 #include <trackbase/InttDefs.h>
0015 #include <trackbase/MvtxDefs.h>
0016 #include <trackbase/TpcDefs.h>
0017 #include <trackbase/TrkrClusterv1.h>
0018 #include <trackbase/TrkrDefs.h>
0019
0020 #include <g4detectors/PHG4CylinderGeom.h>  // for PHG4CylinderGeom
0021 #include <g4detectors/PHG4CylinderGeomContainer.h>
0022 #include <g4detectors/PHG4TpcGeom.h>
0023 #include <g4detectors/PHG4TpcGeomContainer.h>
0024
0025 #include <mvtx/CylinderGeom_Mvtx.h>
0026 #include <mvtx/CylinderGeom_MvtxHelper.h>
0027
0028 #include <intt/CylinderGeomIntt.h>
0029 #include <intt/CylinderGeomInttHelper.h>
0030
0031 #include <phool/PHCompositeNode.h>
0032 #include <phool/getClass.h>
0033 #include <phool/phool.h>  // for PHWHERE
0034
0035 #include <phparameter/PHParameters.h>
0036 #include <phparameter/PHParametersContainer.h>
0037
0038 #include <TVector3.h>
0039
0040 #include <algorithm>
0041 #include <cassert>
0042 #include <cmath>    // for sqrt, fabs
0043 #include <cstdlib>  // for abs
0044 #include <iostream>
0045 #include <limits>
0046 #include <map>
0047 #include <set>
0048 #include <utility>
0049
0050 SvtxTruthEval::SvtxTruthEval(PHCompositeNode* topNode)
0051   : _basetrutheval(topNode)
0052 {
0053   get_node_pointers(topNode);
0054 }
0055
0056 SvtxTruthEval::~SvtxTruthEval()
0057 {
0058   if (_verbosity > 1)
0059   {
0060     if ((_errors > 0) || (_verbosity > 1))
0061     {
0062       std::cout << "SvtxTruthEval::~SvtxTruthEval() - Error Count: " << _errors << std::endl;
0063     }
0064   }
0065 }
0066
0067 void SvtxTruthEval::next_event(PHCompositeNode* topNode)
0068 {
0069   _cache_all_truth_hits.clear();
0070   _cache_all_truth_hits_g4particle.clear();
0071   _cache_all_truth_clusters_g4particle.clear();
0072   _cache_get_innermost_truth_hit.clear();
0073   _cache_get_outermost_truth_hit.clear();
0074   _cache_get_primary_particle_g4hit.clear();
0075
0076   _basetrutheval.next_event(topNode);
0077
0078   get_node_pointers(topNode);
0079 }
0080
0081 /// \todo this copy may be too expensive to call a lot...
0082 std::set<PHG4Hit*> SvtxTruthEval::all_truth_hits()
0083 {
0084   if (!has_node_pointers())
0085   {
0086     ++_errors;
0087     return std::set<PHG4Hit*>();
0088   }
0089
0090   if (_do_cache)
0091   {
0092     if (!_cache_all_truth_hits.empty())
0093     {
0094       return _cache_all_truth_hits;
0095     }
0096   }
0097
0098   // since the SVTX can be composed of two different trackers this is a
0099   // handy function to spill out all the g4hits from both "detectors"
0100
0101   std::set<PHG4Hit*> truth_hits;
0102
0103   // loop over all the g4hits in the cylinder layers
0104   if (_g4hits_svtx)
0105   {
0106     for (PHG4HitContainer::ConstIterator g4iter = _g4hits_svtx->getHits().first;
0107          g4iter != _g4hits_svtx->getHits().second;
0108          ++g4iter)
0109     {
0110       PHG4Hit* g4hit = g4iter->second;
0111       truth_hits.insert(g4hit);
0112     }
0113   }
0114
0115   // loop over all the g4hits in the ladder layers
0116   if (_g4hits_tracker)
0117   {
0118     for (PHG4HitContainer::ConstIterator g4iter = _g4hits_tracker->getHits().first;
0119          g4iter != _g4hits_tracker->getHits().second;
0120          ++g4iter)
0121     {
0122       PHG4Hit* g4hit = g4iter->second;
0123       truth_hits.insert(g4hit);
0124     }
0125   }
0126
0127   // loop over all the g4hits in the maps layers
0128   if (_g4hits_maps)
0129   {
0130     for (PHG4HitContainer::ConstIterator g4iter = _g4hits_maps->getHits().first;
0131          g4iter != _g4hits_maps->getHits().second;
0132          ++g4iter)
0133     {
0134       PHG4Hit* g4hit = g4iter->second;
0135       truth_hits.insert(g4hit);
0136     }
0137   }
0138
0139   // loop over all the g4hits in the maicromegas layers
0140   if (_g4hits_mms)
0141   {
0142     for (PHG4HitContainer::ConstIterator g4iter = _g4hits_mms->getHits().first;
0143          g4iter != _g4hits_mms->getHits().second;
0144          ++g4iter)
0145     {
0146       PHG4Hit* g4hit = g4iter->second;
0147       truth_hits.insert(g4hit);
0148     }
0149   }
0150
0151   if (_do_cache)
0152   {
0153     _cache_all_truth_hits = truth_hits;
0154   }
0155
0156   return truth_hits;
0157 }
0158
0159 std::set<PHG4Hit*> SvtxTruthEval::all_truth_hits(PHG4Particle* particle)
0160 {
0161   if (!has_node_pointers())
0162   {
0163     ++_errors;
0164     return std::set<PHG4Hit*>();
0165   }
0166
0167   if (_strict)
0168   {
0169     assert(particle);
0170   }
0171   else if (!particle)
0172   {
0173     ++_errors;
0174     return std::set<PHG4Hit*>();
0175   }
0176   //  if( _cache_all_truth_hits_g4particle.count(particle)==0){
0177   if (_cache_all_truth_hits_g4particle.empty())
0178   {
0179     FillTruthHitsFromParticleCache();
0180   }
0181
0182   if (_do_cache)
0183   {
0184     std::map<PHG4Particle*, std::set<PHG4Hit*>>::iterator iter =
0185         _cache_all_truth_hits_g4particle.find(particle);
0186     if (iter != _cache_all_truth_hits_g4particle.end())
0187     {
0188       return iter->second;
0189     }
0190   }
0191   // return empty if we dont find anything int he cache
0192
0193   std::set<PHG4Hit*> truth_hits;
0194   return truth_hits;
0195 }
0196
0197 void SvtxTruthEval::FillTruthHitsFromParticleCache()
0198 {
0199   std::multimap<const int, PHG4Hit*> temp_clusters_from_particles;
0200   // loop over all the g4hits in the cylinder layers
0201   if (_g4hits_svtx)
0202   {
0203     for (PHG4HitContainer::ConstIterator g4iter = _g4hits_svtx->getHits().first;
0204          g4iter != _g4hits_svtx->getHits().second;
0205          ++g4iter)
0206     {
0207       PHG4Hit* g4hit = g4iter->second;
0208       temp_clusters_from_particles.insert(std::make_pair(g4hit->get_trkid(), g4hit));
0209     }
0210   }
0211
0212   // loop over all the g4hits in the ladder layers
0213   if (_g4hits_tracker)
0214   {
0215     for (PHG4HitContainer::ConstIterator g4iter = _g4hits_tracker->getHits().first;
0216          g4iter != _g4hits_tracker->getHits().second;
0217          ++g4iter)
0218     {
0219       PHG4Hit* g4hit = g4iter->second;
0220       temp_clusters_from_particles.insert(std::make_pair(g4hit->get_trkid(), g4hit));
0221     }
0222   }
0223
0224   // loop over all the g4hits in the maps ladder layers
0225   if (_g4hits_maps)
0226   {
0227     for (PHG4HitContainer::ConstIterator g4iter = _g4hits_maps->getHits().first;
0228          g4iter != _g4hits_maps->getHits().second;
0229          ++g4iter)
0230     {
0231       PHG4Hit* g4hit = g4iter->second;
0232       temp_clusters_from_particles.insert(std::make_pair(g4hit->get_trkid(), g4hit));
0233     }
0234   }
0235
0236   // loop over all the g4hits in the micromegas layers
0237   if (_g4hits_mms)
0238   {
0239     for (PHG4HitContainer::ConstIterator g4iter = _g4hits_mms->getHits().first;
0240          g4iter != _g4hits_mms->getHits().second;
0241          ++g4iter)
0242     {
0243       PHG4Hit* g4hit = g4iter->second;
0244       temp_clusters_from_particles.insert(std::make_pair(g4hit->get_trkid(), g4hit));
0245     }
0246   }
0247
0248   PHG4TruthInfoContainer::ConstRange range = _truthinfo->GetParticleRange();
0249   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0250        iter != range.second; ++iter)
0251   {
0252     PHG4Particle* g4particle = iter->second;
0253     std::set<PHG4Hit*> truth_hits;
0254     std::multimap<const int, PHG4Hit*>::const_iterator lower_bound = temp_clusters_from_particles.lower_bound(g4particle->get_track_id());
0255     std::multimap<const int, PHG4Hit*>::const_iterator upper_bound = temp_clusters_from_particles.upper_bound(g4particle->get_track_id());
0256     std::multimap<const int, PHG4Hit*>::const_iterator cfp_iter;
0257     for (cfp_iter = lower_bound; cfp_iter != upper_bound; ++cfp_iter)
0258     {
0259       PHG4Hit* g4hit = cfp_iter->second;
0260       truth_hits.insert(g4hit);
0261     }
0262     _cache_all_truth_hits_g4particle.insert(std::make_pair(g4particle, truth_hits));
0263   }
0264 }
0265
0266 std::map<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster>> SvtxTruthEval::all_truth_clusters(PHG4Particle* particle)
0267 {
0268   using output_type_t = std::map<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster>>;
0269   if (!has_node_pointers())
0270   {
0271     ++_errors;
0272     return output_type_t();
0273   }
0274
0275   if (_strict)
0276   {
0277     assert(particle);
0278   }
0279   else if (!particle)
0280   {
0281     ++_errors;
0282     return output_type_t();
0283   }
0284
0285   if (_do_cache)
0286   {
0287     const auto iter = _cache_all_truth_clusters_g4particle.find(particle);
0288     if (iter != _cache_all_truth_clusters_g4particle.end())
0289     {
0290       return iter->second;
0291     }
0292   }
0293
0294   if (_verbosity > 1)
0295   {
0296     std::cout << PHWHERE << " Truth clustering for particle " << particle->get_track_id() << std::endl;
0297   };
0298
0299   // get all g4hits for this particle
0300   std::set<PHG4Hit*> g4hits = all_truth_hits(particle);
0301
0302   float ng4hits = g4hits.size();
0303   if (ng4hits == 0)
0304   {
0305     return output_type_t();
0306   }
0307
0308   // container for storing truth clusters  //std::map<unsigned int, TrkrCluster*> truth_clusters;
0309   output_type_t truth_clusters;
0310
0311   // convert truth hits for this particle to truth clusters in each layer
0312   // loop over layers
0313
0314   unsigned int layer;
0315   for (layer = 0; layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms; ++layer)
0316   {
0317     float gx = std::numeric_limits<float>::quiet_NaN();
0318     float gy = std::numeric_limits<float>::quiet_NaN();
0319     float gz = std::numeric_limits<float>::quiet_NaN();
0320     float gt = std::numeric_limits<float>::quiet_NaN();
0321     float gedep = std::numeric_limits<float>::quiet_NaN();
0322
0323     std::vector<PHG4Hit*> contributing_hits;
0324     std::vector<double> contributing_hits_energy;
0325     std::vector<std::vector<double>> contributing_hits_entry;
0326     std::vector<std::vector<double>> contributing_hits_exit;
0327     LayerClusterG4Hits(g4hits, contributing_hits, contributing_hits_energy, contributing_hits_entry, contributing_hits_exit, layer, gx, gy, gz, gt, gedep);
0328     if (!(gedep > 0))
0329     {
0330       continue;
0331     }
0332
0333     // we have the cluster in this layer from this truth particle
0334     // add the cluster to a TrkrCluster object
0335     TrkrDefs::cluskey ckey;
0336     if (layer >= _nlayers_maps + _nlayers_intt && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)  // in TPC
0337     {
0338       unsigned int side = 0;
0339       if (gz > 0)
0340       {
0341         side = 1;
0342       }
0343       // need dummy sector here
0344       unsigned int sector = 0;
0345       ckey = TpcDefs::genClusKey(layer, sector, side, iclus);
0346     }
0347     else if (layer < _nlayers_maps)  // in MVTX
0348     {
0349       unsigned int stave = 0;
0350       unsigned int chip = 0;
0351       unsigned int strobe = 0;
0352       ckey = MvtxDefs::genClusKey(layer, stave, chip, strobe, iclus);
0353     }
0354     else if (layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt)  // in INTT
0355     {
0356       // dummy ladder and phi ID
0357       unsigned int ladderzid = 0;
0358       unsigned int ladderphiid = 0;
0359       uint16_t crossing = 0;
0360       ckey = InttDefs::genClusKey(layer, ladderzid, ladderphiid, crossing, iclus);
0361     }
0362     else if (layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)  // in MICROMEGAS
0363     {
0364       unsigned int tile = 0;
0365       MicromegasDefs::SegmentationType segtype;
0366       segtype = MicromegasDefs::SegmentationType::SEGMENTATION_PHI;
0367       TrkrDefs::hitsetkey hkey = MicromegasDefs::genHitSetKey(layer, segtype, tile);
0368       ckey = TrkrDefs::genClusKey(hkey, iclus);
0369     }
0370     else
0371     {
0372       std::cout << PHWHERE << "Bad layer number: " << layer << std::endl;
0373       continue;
0374     }
0375
0376     auto clus = std::make_shared<TrkrClusterv1>();
0377     iclus++;
0378
0379     // estimate cluster ADC value
0380     unsigned int adc_value = getAdcValue(gedep);
0381     // std::cout << " gedep " << gedep << " adc_value " << adc_value << std::endl;
0382     clus->setAdc(adc_value);
0383     clus->setPosition(0, gx);
0384     clus->setPosition(1, gy);
0385     clus->setPosition(2, gz);
0386     clus->setGlobal();
0387
0388     // record which g4hits contribute to this truth cluster
0389     for (auto& contributing_hit : contributing_hits)
0390     {
0391       _truth_cluster_truth_hit_map.insert(std::make_pair(ckey, contributing_hit));
0392     }
0393
0394     // Estimate the size of the truth cluster
0395     float g4phisize = std::numeric_limits<float>::quiet_NaN();
0396     float g4zsize = std::numeric_limits<float>::quiet_NaN();
0397     G4ClusterSize(ckey, layer, contributing_hits_entry, contributing_hits_exit, g4phisize, g4zsize);
0398
0399     for (int i1 = 0; i1 < 3; ++i1)
0400     {
0401       for (int i2 = 0; i2 < 3; ++i2)
0402       {
0403         clus->setSize(i1, i2, 0.0);
0404         clus->setError(i1, i2, 0.0);
0405       }
0406     }
0407     clus->setError(0, 0, gedep);  // stores truth energy
0408     clus->setSize(1, 1, g4phisize);
0409     clus->setSize(2, 2, g4zsize);
0410     clus->setError(1, 1, g4phisize / std::sqrt(12));
0411     clus->setError(2, 2, g4zsize / std::sqrt(12.0));
0412
0413     truth_clusters.insert(std::make_pair(ckey, clus));
0414
0415   }  // end loop over layers for this particle
0416
0417   if (_do_cache)
0418   {
0419     _cache_all_truth_clusters_g4particle.insert(std::make_pair(particle, truth_clusters));
0420   }
0421
0422   return truth_clusters;
0423 }
0424
0425 void SvtxTruthEval::LayerClusterG4Hits(const std::set<PHG4Hit*>& truth_hits, std::vector<PHG4Hit*>& contributing_hits, std::vector<double>& contributing_hits_energy, std::vector<std::vector<double>>& contributing_hits_entry, std::vector<std::vector<double>>& contributing_hits_exit, float layer, float& x, float& y, float& z, float& t, float& e)
0426 {
0427   bool use_geo = true;
0428
0429   // Given a set of g4hits, cluster them within a given layer of the TPC
0430
0431   float gx = 0.0;
0432   float gy = 0.0;
0433   float gz = 0.0;
0434   float gr = 0.0;
0435   float gt = 0.0;
0436   float gwt = 0.0;
0437
0438   if (layer >= _nlayers_maps + _nlayers_intt && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)  // in TPC
0439   {
0440     // std::cout << "layer = " << layer << " _nlayers_maps " << _nlayers_maps << " _nlayers_intt " << _nlayers_intt << std::endl;
0441
0442     // This calculates the truth cluster position for the TPC from all of the contributing g4hits from a g4particle, typically 2-4 for the TPC
0443     // Complicated, since only the part of the energy that is collected within a layer contributes to the position
0444     //===============================================================================
0445
0446     PHG4TpcGeom* GeoLayer = _tpc_geom_container->GetLayerCellGeom(layer);
0447     // get layer boundaries here for later use
0448     // radii of layer boundaries
0449     float rbin = GeoLayer->get_radius() - GeoLayer->get_thickness() / 2.0;
0450     float rbout = GeoLayer->get_radius() + GeoLayer->get_thickness() / 2.0;
0451
0452     if (_verbosity > 1)
0453     {
0454       std::cout << " TruthEval::LayerCluster hits for layer  " << layer << " with rbin " << rbin << " rbout " << rbout << std::endl;
0455     }
0456
0457     // we do not assume that the truth hits know what layer they are in
0458     for (auto *this_g4hit : truth_hits)
0459     {
0460       float rbegin = std::sqrt(this_g4hit->get_x(0) * this_g4hit->get_x(0) + this_g4hit->get_y(0) * this_g4hit->get_y(0));
0461       float rend = std::sqrt(this_g4hit->get_x(1) * this_g4hit->get_x(1) + this_g4hit->get_y(1) * this_g4hit->get_y(1));
0462       // std::cout << " Eval: g4hit " << this_g4hit->get_hit_id() <<  " layer " << layer << " rbegin " << rbegin << " rend " << rend << std::endl;
0463
0464       // make sure the entry point is at lower radius
0465       float xl[2];
0466       float yl[2];
0467       float zl[2];
0468
0469       if (rbegin < rend)
0470       {
0471         xl[0] = this_g4hit->get_x(0);
0472         yl[0] = this_g4hit->get_y(0);
0473         zl[0] = this_g4hit->get_z(0);
0474         xl[1] = this_g4hit->get_x(1);
0475         yl[1] = this_g4hit->get_y(1);
0476         zl[1] = this_g4hit->get_z(1);
0477       }
0478       else
0479       {
0480         xl[0] = this_g4hit->get_x(1);
0481         yl[0] = this_g4hit->get_y(1);
0482         zl[0] = this_g4hit->get_z(1);
0483         xl[1] = this_g4hit->get_x(0);
0484         yl[1] = this_g4hit->get_y(0);
0485         zl[1] = this_g4hit->get_z(0);
0486         std::swap(rbegin, rend);
0487         // std::cout << "swapped in and out " << std::endl;
0488       }
0489
0490       // check that the g4hit is not completely outside the cluster layer. Just skip this g4hit if it is
0491       if (rbegin < rbin && rend < rbin)
0492       {
0493         continue;
0494       }
0495       if (rbegin > rbout && rend > rbout)
0496       {
0497         continue;
0498       }
0499
0500       if (_verbosity > 1)
0501       {
0502         std::cout << "     keep g4hit with rbegin " << rbegin << " rend " << rend
0503                   << "         xbegin " << xl[0] << " xend " << xl[1]
0504                   << " ybegin " << yl[0] << " yend " << yl[1]
0505                   << " zbegin " << zl[0] << " zend " << zl[1]
0506                   << std::endl;
0507       }
0508
0509       float xin = xl[0];
0510       float yin = yl[0];
0511       float zin = zl[0];
0512       float xout = xl[1];
0513       float yout = yl[1];
0514       float zout = zl[1];
0515
0516       float local_t = std::numeric_limits<float>::quiet_NaN();
0517
0518       if (rbegin < rbin)
0519       {
0520         // line segment begins before boundary, find where it crosses
0521         local_t = line_circle_intersection(xl, yl, zl, rbin);
0522         if (local_t > 0)
0523         {
0524           xin = xl[0] + local_t * (xl[1] - xl[0]);
0525           yin = yl[0] + local_t * (yl[1] - yl[0]);
0526           zin = zl[0] + local_t * (zl[1] - zl[0]);
0527         }
0528       }
0529
0530       if (rend > rbout)
0531       {
0532         // line segment ends after boundary, find where it crosses
0533         local_t = line_circle_intersection(xl, yl, zl, rbout);
0534         if (local_t > 0)
0535         {
0536           xout = xl[0] + local_t * (xl[1] - xl[0]);
0537           yout = yl[0] + local_t * (yl[1] - yl[0]);
0538           zout = zl[0] + local_t * (zl[1] - zl[0]);
0539         }
0540       }
0541
0542       double rin = std::sqrt(xin * xin + yin * yin);
0543       double rout = std::sqrt(xout * xout + yout * yout);
0544
0545       // we want only the fraction of edep inside the layer
0546       double efrac = this_g4hit->get_edep() * (rout - rin) / (rend - rbegin);
0547       gx += (xin + xout) * 0.5 * efrac;
0548       gy += (yin + yout) * 0.5 * efrac;
0549       gz += (zin + zout) * 0.5 * efrac;
0550       gt += this_g4hit->get_avg_t() * efrac;
0551       gr += (rin + rout) * 0.5 * efrac;
0552       gwt += efrac;
0553
0554       if (_verbosity > 1)
0555       {
0556         std::cout << "      rin  " << rin << " rout " << rout
0557                   << " xin " << xin << " xout " << xout << " yin " << yin << " yout " << yout << " zin " << zin << " zout " << zout
0558                   << " edep " << this_g4hit->get_edep()
0559                   << " this_edep " << efrac << std::endl;
0560         std::cout << "              xavge " << (xin + xout) * 0.5 << " yavge " << (yin + yout) * 0.5 << " zavge " << (zin + zout) * 0.5 << " ravge " << (rin + rout) * 0.5
0561                   << std::endl;
0562       }
0563
0564       // Capture entry and exit points
0565       std::vector<double> entry_loc;
0566       entry_loc.push_back(xin);
0567       entry_loc.push_back(yin);
0568       entry_loc.push_back(zin);
0569       std::vector<double> exit_loc;
0570       exit_loc.push_back(xout);
0571       exit_loc.push_back(yout);
0572       exit_loc.push_back(zout);
0573
0574       // this_g4hit is inside the layer, add it to the vectors
0575       contributing_hits.push_back(this_g4hit);
0576       contributing_hits_energy.push_back(this_g4hit->get_edep() * (zout - zin) / (zl[1] - zl[0]));
0577       contributing_hits_entry.push_back(entry_loc);
0578       contributing_hits_exit.push_back(exit_loc);
0579
0580     }  // loop over contributing hits
0581
0582     if (gwt == 0)
0583     {
0584       e = gwt;
0585       return;  // will be discarded
0586     }
0587
0588     gx /= gwt;
0589     gy /= gwt;
0590     gz /= gwt;
0591     gr = (rbin + rbout) * 0.5;
0592     gt /= gwt;
0593
0594     if (_verbosity > 1)
0595     {
0596       std::cout << " weighted means:   gx " << gx << " gy " << gy << " gz " << gz << " gr " << gr << " e " << gwt << std::endl;
0597     }
0598
0599     if (use_geo)
0600     {
0601       // The energy weighted values above have significant scatter due to fluctuations in the energy deposit from Geant
0602       // Calculate the geometric mean positions instead
0603       float rentry = 999.0;
0604       float xentry = 999.0;
0605       float yentry = 999.0;
0606       float zentry = 999.0;
0607       float rexit = -999.0;
0608       float xexit = -999.0;
0609       float yexit = -999.0;
0610       float zexit = -999.0;
0611
0612       for (unsigned int ientry = 0; ientry < contributing_hits_entry.size(); ++ientry)
0613       {
0614         float tmpx = contributing_hits_entry[ientry][0];
0615         float tmpy = contributing_hits_entry[ientry][1];
0616         float tmpr = std::sqrt(tmpx * tmpx + tmpy * tmpy);
0617
0618         if (tmpr < rentry)
0619         {
0620           rentry = tmpr;
0621           xentry = contributing_hits_entry[ientry][0];
0622           yentry = contributing_hits_entry[ientry][1];
0623           zentry = contributing_hits_entry[ientry][2];
0624         }
0625
0626         tmpx = contributing_hits_exit[ientry][0];
0627         tmpy = contributing_hits_exit[ientry][1];
0628         tmpr = std::sqrt(tmpx * tmpx + tmpy * tmpy);
0629
0630         if (tmpr > rexit)
0631         {
0632           rexit = tmpr;
0633           xexit = contributing_hits_exit[ientry][0];
0634           yexit = contributing_hits_exit[ientry][1];
0635           zexit = contributing_hits_exit[ientry][2];
0636         }
0637       }
0638
0639       float geo_r = (rentry + rexit) * 0.5;
0640       float geo_x = (xentry + xexit) * 0.5;
0641       float geo_y = (yentry + yexit) * 0.5;
0642       float geo_z = (zentry + zexit) * 0.5;
0643
0644       if (rexit > 0)
0645       {
0646         gx = geo_x;
0647         gy = geo_y;
0648         gz = geo_z;
0649         gr = geo_r;
0650       }
0651
0652       if (_verbosity > 1)
0653       {
0654         std::cout << "      rentry  " << rentry << " rexit " << rexit
0655                   << " xentry " << xentry << " xexit " << xexit << " yentry " << yentry << " yexit " << yexit << " zentry " << zentry << " zexit " << zexit << std::endl;
0656
0657         std::cout << " geometric means: geo_x " << geo_x << " geo_y " << geo_y << " geo_z " << geo_z << " geo r " << geo_r << " e " << gwt << std::endl
0658                   << std::endl;
0659       }
0660     }
0661
0662   }  // if TPC
0663   else
0664   {
0665     // not TPC, one g4hit per cluster
0666     for (auto *this_g4hit : truth_hits)
0667     {
0668       if (this_g4hit->get_layer() != (unsigned int) layer)
0669       {
0670         continue;
0671       }
0672
0673       gx = this_g4hit->get_avg_x();
0674       gy = this_g4hit->get_avg_y();
0675       gz = this_g4hit->get_avg_z();
0676       gt = this_g4hit->get_avg_t();
0677       gwt += this_g4hit->get_edep();
0678
0679       // Capture entry and exit points
0680       std::vector<double> entry_loc;
0681       entry_loc.push_back(this_g4hit->get_x(0));
0682       entry_loc.push_back(this_g4hit->get_y(0));
0683       entry_loc.push_back(this_g4hit->get_z(0));
0684       std::vector<double> exit_loc;
0685       exit_loc.push_back(this_g4hit->get_x(1));
0686       exit_loc.push_back(this_g4hit->get_y(1));
0687       exit_loc.push_back(this_g4hit->get_z(1));
0688
0689       // this_g4hit is inside the layer, add it to the vectors
0690       contributing_hits.push_back(this_g4hit);
0691       contributing_hits_energy.push_back(this_g4hit->get_edep());
0692       contributing_hits_entry.push_back(entry_loc);
0693       contributing_hits_exit.push_back(exit_loc);
0694     }
0695   }  // not TPC
0696
0697   x = gx;
0698   y = gy;
0699   z = gz;
0700   t = gt;
0701   e = gwt;
0702
0703   return;
0704 }
0705
0706 void SvtxTruthEval::G4ClusterSize(TrkrDefs::cluskey ckey, unsigned int layer, const std::vector<std::vector<double>> &contributing_hits_entry, const std::vector<std::vector<double>> &contributing_hits_exit, float& g4phisize, float& g4zsize)
0707 {
0708   // sort the contributing g4hits in radius
0709   double inner_radius = 100.;
0710   double inner_x = std::numeric_limits<double>::quiet_NaN();
0711   double inner_y = std::numeric_limits<double>::quiet_NaN();
0712   double inner_z = std::numeric_limits<double>::quiet_NaN();
0713
0714   double outer_radius = 0.;
0715   double outer_x = std::numeric_limits<double>::quiet_NaN();
0716   double outer_y = std::numeric_limits<double>::quiet_NaN();
0717   double outer_z = std::numeric_limits<double>::quiet_NaN();
0718
0719   for (unsigned int ihit = 0; ihit < contributing_hits_entry.size(); ++ihit)
0720   {
0721     double rad1 = std::sqrt(pow(contributing_hits_entry[ihit][0], 2) + pow(contributing_hits_entry[ihit][1], 2));
0722     if (rad1 < inner_radius)
0723     {
0724       inner_radius = rad1;
0725       inner_x = contributing_hits_entry[ihit][0];
0726       inner_y = contributing_hits_entry[ihit][1];
0727       inner_z = contributing_hits_entry[ihit][2];
0728     }
0729
0730     double rad2 = std::sqrt(pow(contributing_hits_exit[ihit][0], 2) + pow(contributing_hits_exit[ihit][1], 2));
0731     if (rad2 > outer_radius)
0732     {
0733       outer_radius = rad2;
0734       outer_x = contributing_hits_exit[ihit][0];
0735       outer_y = contributing_hits_exit[ihit][1];
0736       outer_z = contributing_hits_exit[ihit][2];
0737     }
0738   }
0739
0740   double inner_phi = atan2(inner_y, inner_x);
0741   double outer_phi = atan2(outer_y, outer_x);
0742   double avge_z = (outer_z + inner_z) / 2.0;
0743
0744   unsigned int side = 0;
0745   if (avge_z < 0)
0746   {
0747     side = 1;
0748   }
0749
0750   // Now fold these with the expected diffusion and shaping widths
0751   // assume spread is +/- equals this many sigmas times diffusion and shaping when extending the size
0752   double sigmas = 2.0;
0753
0754   double radius = (inner_radius + outer_radius) / 2.;
0755   if (radius > 28 && radius < 80)  // TPC
0756   {
0757     PHG4TpcGeom* layergeom = _tpc_geom_container->GetLayerCellGeom(layer);
0758
0759     double tpc_max_driftlength = layergeom->get_max_driftlength();
0760     double drift_velocity = layergeom->get_drift_velocity_sim();
0761
0762     // Phi size
0763     //======
0764     double diffusion_trans = 0.006;  // cm/SQRT(cm)
0765     double phidiffusion = diffusion_trans * std::sqrt(tpc_max_driftlength - fabs(avge_z));
0766
0767     double added_smear_trans = 0.085;  // cm
0768     double gem_spread = 0.04;          // 400 microns
0769
0770     if (outer_phi < inner_phi)
0771     {
0772       std::swap(outer_phi, inner_phi);
0773     }
0774
0775     // convert diffusion from cm to radians
0776     double g4max_phi = outer_phi + sigmas * std::sqrt(pow(phidiffusion, 2) + pow(added_smear_trans, 2) + pow(gem_spread, 2)) / radius;
0777     double g4min_phi = inner_phi - sigmas * std::sqrt(pow(phidiffusion, 2) + pow(added_smear_trans, 2) + pow(gem_spread, 2)) / radius;
0778
0779     unsigned int phibinwidth = 1;
0780     if(std::isnan(g4min_phi) || std::isnan(g4max_phi))
0781       {
0782       if (_verbosity > 1)
0783     {
0784       std::cout << " g4min_phi " << g4min_phi << " g4max_phi " << g4max_phi << std::endl;
0785       std::cout << "     outer_phi " << outer_phi << " inner_phi " << inner_phi << " radius " << radius << std::endl;
0786     }
0787       }
0788     else
0789       {
0790     // find the bins containing these max and min z edges
0791     unsigned int phibinmin = layergeom->get_phibin(g4min_phi, side);
0792     unsigned int phibinmax = layergeom->get_phibin(g4max_phi, side);
0793     phibinwidth = phibinmax - phibinmin + 1;
0794       }
0795     g4phisize = (double) phibinwidth * layergeom->get_phistep() * layergeom->get_radius();
0796
0797     // Z size
0798     //=====
0799     double g4max_z = 0;
0800     double g4min_z = 0;
0801
0802     outer_z = fabs(outer_z);
0803     inner_z = fabs(inner_z);
0804
0805     double diffusion_long = 0.015;  // cm/SQRT(cm)
0806     double zdiffusion = diffusion_long * std::sqrt(tpc_max_driftlength - fabs(avge_z));
0807     double zshaping_lead = 32.0 * drift_velocity;  // ns * cm/ns = cm
0808     double zshaping_tail = 48.0 * drift_velocity;
0809     double added_smear_long = 0.105;  // cm
0810
0811     // largest z reaches gems first, make that the outer z
0812     if (outer_z < inner_z)
0813     {
0814       std::swap(outer_z, inner_z);
0815     }
0816     g4max_z = outer_z + sigmas * std::sqrt(pow(zdiffusion, 2) + pow(added_smear_long, 2) + pow(zshaping_lead, 2));
0817     g4min_z = inner_z - sigmas * std::sqrt(pow(zdiffusion, 2) + pow(added_smear_long, 2) + pow(zshaping_tail, 2));
0818
0819     // find the bins containing these max and min z edges
0820     unsigned int binmin = layergeom->get_zbin(g4min_z);
0821     unsigned int binmax = layergeom->get_zbin(g4max_z);
0822     if (binmax < binmin)
0823     {
0824       std::swap(binmax, binmin);
0825     }
0826     unsigned int binwidth = binmax - binmin + 1;
0827
0828     // multiply total number of bins that include the edges by the bin size
0829     g4zsize = (double) binwidth * layergeom->get_zstep();
0830   }
0831   else if (radius > 6 && radius < 20)  // INTT
0832   {
0833     // All we have is the position and layer number
0834
0835     CylinderGeomIntt* layergeom = dynamic_cast<CylinderGeomIntt*>(_intt_geom_container->GetLayerGeom(layer));
0836
0837     // inner location
0838     double world_inner[3] = {inner_x, inner_y, inner_z};
0839     TVector3 world_inner_vec = {inner_x, inner_y, inner_z};
0840
0841     int segment_z_bin;
0842     int segment_phi_bin;
0843     layergeom->find_indices_from_world_location(segment_z_bin, segment_phi_bin, world_inner);
0844
0845     TrkrDefs::hitsetkey hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(ckey);
0846     auto surf = _tgeometry->maps().getSiliconSurface(hitsetkey);
0847     TVector3 local_inner_vec = CylinderGeomInttHelper::get_local_from_world_coords(surf, _tgeometry, world_inner);
0848     double yin = local_inner_vec[1];
0849     double zin = local_inner_vec[2];
0850     int strip_y_index;
0851     int strip_z_index;
0852     layergeom->find_strip_index_values(segment_z_bin, yin, zin, strip_y_index, strip_z_index);
0853
0854     // outer location
0855     double world_outer[3] = {outer_x, outer_y, outer_z};
0856     TVector3 world_outer_vec = {outer_x, outer_y, outer_z};
0857
0858     layergeom->find_indices_from_world_location(segment_z_bin, segment_phi_bin, world_outer);
0859     TrkrDefs::hitsetkey ohitsetkey = TrkrDefs::getHitSetKeyFromClusKey(ckey);
0860     auto osurf = _tgeometry->maps().getSiliconSurface(ohitsetkey);
0861     TVector3 local_outer_vec = CylinderGeomInttHelper::get_local_from_world_coords(osurf, _tgeometry, world_outer_vec);
0862     double yout = local_outer_vec[1];
0863     double zout = local_outer_vec[2];
0864     int strip_y_index_out;
0865     int strip_z_index_out;
0866     layergeom->find_strip_index_values(segment_z_bin, yout, zout, strip_y_index_out, strip_z_index_out);
0867
0868     int strips = abs(strip_y_index_out - strip_y_index) + 1;
0869     int cols = abs(strip_z_index_out - strip_z_index) + 1;
0870
0871     double strip_width = (double) strips * layergeom->get_strip_y_spacing();  // cm
0872     double strip_length = (double) cols * layergeom->get_strip_z_spacing();   // cm
0873
0874     g4phisize = strip_width;
0875     g4zsize = strip_length;
0876
0877     /*
0878     if(Verbosity() > 1)
0879       std::cout << " INTT: layer " << layer << " strips " << strips << " strip pitch " <<  layergeom->get_strip_y_spacing() << " g4phisize "<< g4phisize
0880            << " columns " << cols << " strip_z_spacing " <<  layergeom->get_strip_z_spacing() << " g4zsize " << g4zsize << std::endl;
0881     */
0882   }
0883   else if (radius > 80)  // MICROMEGAS
0884   {
0885     // made up for now
0886     g4phisize = 300e-04;
0887     g4zsize = 300e-04;
0888   }
0889   else  // MVTX
0890   {
0891     unsigned int stave;
0892     unsigned int stave_outer;
0893     unsigned int chip;
0894     unsigned int chip_outer;
0895     int row;
0896     int row_outer;
0897     int column;
0898     int column_outer;
0899
0900     // add diffusion to entry and exit locations
0901     double max_diffusion_radius = 25.0e-4;  // 25 microns
0902     double min_diffusion_radius = 8.0e-4;   // 8 microns
0903
0904     CylinderGeom_Mvtx* layergeom = dynamic_cast<CylinderGeom_Mvtx*>(_mvtx_geom_container->GetLayerGeom(layer));
0905
0906     TVector3 world_inner = {inner_x, inner_y, inner_z};
0907     std::vector<double> world_inner_vec = {world_inner[0], world_inner[1], world_inner[2]};
0908     layergeom->get_sensor_indices_from_world_coords(world_inner_vec, stave, chip);
0909     TrkrDefs::hitsetkey ihitsetkey = TrkrDefs::getHitSetKeyFromClusKey(ckey);
0910     auto isurf = _tgeometry->maps().getSiliconSurface(ihitsetkey);
0911     TVector3 local_inner = CylinderGeom_MvtxHelper::get_local_from_world_coords(isurf, _tgeometry, world_inner);
0912
0913     TVector3 world_outer = {outer_x, outer_y, outer_z};
0914     std::vector<double> world_outer_vec = {world_outer[0], world_outer[1], world_outer[2]};
0915     layergeom->get_sensor_indices_from_world_coords(world_outer_vec, stave_outer, chip_outer);
0916     TrkrDefs::hitsetkey ohitsetkey = TrkrDefs::getHitSetKeyFromClusKey(ckey);
0917     auto osurf = _tgeometry->maps().getSiliconSurface(ohitsetkey);
0918     TVector3 local_outer = CylinderGeom_MvtxHelper::get_local_from_world_coords(osurf, _tgeometry, world_outer);
0919
0920     double diff = max_diffusion_radius * 0.6;  // factor of 0.6 gives decent agreement with low occupancy reco clusters
0921     if (local_outer[0] < local_inner[0])
0922     {
0923       diff = -diff;
0924     }
0925     local_outer[0] += diff;
0926     local_inner[0] -= diff;
0927
0928     double diff_outer = min_diffusion_radius * 0.6;
0929     if (local_outer[2] < local_inner[2])
0930     {
0931       diff_outer = -diff_outer;
0932     }
0933     local_outer[2] += diff_outer;
0934     local_inner[2] -= diff_outer;
0935
0936     layergeom->get_pixel_from_local_coords(local_inner, row, column);
0937     layergeom->get_pixel_from_local_coords(local_outer, row_outer, column_outer);
0938
0939     if (row_outer < row)
0940     {
0941       std::swap(row_outer, row);
0942     }
0943     unsigned int rows = row_outer - row + 1;
0944     g4phisize = (double) rows * layergeom->get_pixel_x();
0945
0946     if (column_outer < column)
0947     {
0948       std::swap(column_outer, column);
0949     }
0950     unsigned int columns = column_outer - column + 1;
0951     g4zsize = (double) columns * layergeom->get_pixel_z();
0952
0953     /*
0954     if(Verbosity() > 1)
0955       std::cout << " MVTX: layer " << layer << " rows " << rows << " pixel x " <<  layergeom->get_pixel_x() << " g4phisize "<< g4phisize
0956            << " columns " << columns << " pixel_z " <<  layergeom->get_pixel_z() << " g4zsize " << g4zsize << std::endl;
0957     */
0958   }
0959 }
0960
0961 std::set<PHG4Hit*> SvtxTruthEval::get_truth_hits_from_truth_cluster(TrkrDefs::cluskey ckey)
0962 {
0963   std::set<PHG4Hit*> g4hit_set;
0964
0965   std::pair<std::multimap<TrkrDefs::cluskey, PHG4Hit*>::iterator,
0966             std::multimap<TrkrDefs::cluskey, PHG4Hit*>::iterator>
0967       ret = _truth_cluster_truth_hit_map.equal_range(ckey);
0968
0969   for (std::multimap<TrkrDefs::cluskey, PHG4Hit*>::iterator jter = ret.first; jter != ret.second; ++jter)
0970   {
0971     g4hit_set.insert(jter->second);
0972   }
0973
0974   return g4hit_set;
0975 }
0976
0977 PHG4Hit* SvtxTruthEval::get_innermost_truth_hit(PHG4Particle* particle)
0978 {
0979   if (!has_node_pointers())
0980   {
0981     ++_errors;
0982     return nullptr;
0983   }
0984
0985   if (_strict)
0986   {
0987     assert(particle);
0988   }
0989   else if (!particle)
0990   {
0991     ++_errors;
0992     return nullptr;
0993   }
0994
0995   PHG4Hit* innermost_hit = nullptr;
0996   float innermost_radius = std::numeric_limits<float>::max();
0997
0998   std::set<PHG4Hit*> truth_hits = all_truth_hits(particle);
0999   for (auto *candidate : truth_hits)
1000   {
1001     float x = candidate->get_x(0);  // use entry points
1002     float y = candidate->get_y(0);  // use entry points
1003     float r = std::sqrt(x * x + y * y);
1004     if (r < innermost_radius)
1005     {
1006       innermost_radius = r;
1007       innermost_hit = candidate;
1008     }
1009   }
1010
1011   return innermost_hit;
1012 }
1013
1014 PHG4Hit* SvtxTruthEval::get_outermost_truth_hit(PHG4Particle* particle)
1015 {
1016   if (!has_node_pointers())
1017   {
1018     ++_errors;
1019     return nullptr;
1020   }
1021
1022   if (_strict)
1023   {
1024     assert(particle);
1025   }
1026   else if (!particle)
1027   {
1028     ++_errors;
1029     return nullptr;
1030   }
1031
1032   PHG4Hit* outermost_hit = nullptr;
1033   float outermost_radius = std::numeric_limits<float>::min();
1034
1035   if (_do_cache)
1036   {
1037     std::map<PHG4Particle*, PHG4Hit*>::iterator iter =
1038         _cache_get_outermost_truth_hit.find(particle);
1039     if (iter != _cache_get_outermost_truth_hit.end())
1040     {
1041       return iter->second;
1042     }
1043   }
1044
1045   std::set<PHG4Hit*> truth_hits = all_truth_hits(particle);
1046   for (auto *candidate : truth_hits)
1047   {
1048     float x = candidate->get_x(1);  // use exit points
1049     float y = candidate->get_y(1);  // use exit points
1050     float r = std::sqrt(x * x + y * y);
1051     if (r > outermost_radius)
1052     {
1053       outermost_radius = r;
1054       outermost_hit = candidate;
1055     }
1056   }
1057   if (_do_cache)
1058   {
1059     _cache_get_outermost_truth_hit.insert(std::make_pair(particle, outermost_hit));
1060   }
1061
1062   return outermost_hit;
1063 }
1064
1065 PHG4Particle* SvtxTruthEval::get_particle(PHG4Hit* g4hit)
1066 {
1067   return _basetrutheval.get_particle(g4hit);
1068 }
1069
1070 int SvtxTruthEval::get_embed(PHG4Particle* particle)
1071 {
1072   return _basetrutheval.get_embed(particle);
1073 }
1074
1075 PHG4VtxPoint* SvtxTruthEval::get_vertex(PHG4Particle* particle)
1076 {
1077   return _basetrutheval.get_vertex(particle);
1078 }
1079
1080 bool SvtxTruthEval::is_primary(PHG4Particle* particle)
1081 {
1082   return _basetrutheval.is_primary(particle);
1083 }
1084
1085 PHG4Particle* SvtxTruthEval::get_primary_particle(PHG4Hit* g4hit)
1086 {
1087   if (!has_node_pointers())
1088   {
1089     ++_errors;
1090     return nullptr;
1091   }
1092
1093   if (_strict)
1094   {
1095     assert(g4hit);
1096   }
1097   else if (!g4hit)
1098   {
1099     ++_errors;
1100     return nullptr;
1101   }
1102
1103   if (_do_cache)
1104   {
1105     std::map<PHG4Hit*, PHG4Particle*>::iterator iter =
1106         _cache_get_primary_particle_g4hit.find(g4hit);
1107     if (iter != _cache_get_primary_particle_g4hit.end())
1108     {
1109       return iter->second;
1110     }
1111   }
1112
1113   PHG4Particle* primary = _basetrutheval.get_primary_particle(g4hit);
1114
1115   if (_do_cache)
1116   {
1117     _cache_get_primary_particle_g4hit.insert(std::make_pair(g4hit, primary));
1118   }
1119
1120   if (_strict)
1121   {
1122     assert(primary);
1123   }
1124   else if (!primary)
1125   {
1126     ++_errors;
1127   }
1128
1129   return primary;
1130 }
1131
1132 PHG4Particle* SvtxTruthEval::get_primary_particle(PHG4Particle* particle)
1133 {
1134   return _basetrutheval.get_primary_particle(particle);
1135 }
1136
1137 PHG4Particle* SvtxTruthEval::get_particle(const int trackid)
1138 {
1139   return _basetrutheval.get_particle(trackid);
1140 }
1141
1142 bool SvtxTruthEval::is_g4hit_from_particle(PHG4Hit* g4hit, PHG4Particle* particle)
1143 {
1144   return _basetrutheval.is_g4hit_from_particle(g4hit, particle);
1145 }
1146
1147 bool SvtxTruthEval::are_same_particle(PHG4Particle* p1, PHG4Particle* p2)
1148 {
1149   return _basetrutheval.are_same_particle(p1, p2);
1150 }
1151
1152 bool SvtxTruthEval::are_same_vertex(PHG4VtxPoint* vtx1, PHG4VtxPoint* vtx2)
1153 {
1154   return _basetrutheval.are_same_vertex(vtx1, vtx2);
1155 }
1156
1157 void SvtxTruthEval::get_node_pointers(PHCompositeNode* topNode)
1158 {
1159   _tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
1160   _truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1161
1162   _g4hits_mms = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
1163   _g4hits_svtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
1164   _g4hits_tracker = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
1165   _g4hits_maps = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
1166
1167   _mms_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
1168   _tpc_geom_container = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
1169   _intt_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
1170   _mvtx_geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
1171
1172   return;
1173 }
1174
1175 bool SvtxTruthEval::has_node_pointers()
1176 {
1177   if (_strict)
1178   {
1179     assert(_truthinfo);
1180   }
1181   else if (!_truthinfo)
1182   {
1183     return false;
1184   }
1185
1186   if (_strict)
1187   {
1188     assert(_g4hits_mms || _g4hits_svtx || _g4hits_tracker || _g4hits_maps);
1189   }
1190   else if (!_g4hits_mms && !_g4hits_svtx && !_g4hits_tracker && !_g4hits_maps)
1191   {
1192     return false;
1193   }
1194
1195   return true;
1196 }
1197
1198 float SvtxTruthEval::line_circle_intersection(float x[], float y[], float z[], float radius)
1199 {
1200   // parameterize the line in terms of t (distance along the line segment, from 0-1) as
1201   // x = x0 + t * (x1-x0); y=y0 + t * (y1-y0); z = z0 + t * (z1-z0)
1202   // parameterize the cylinder (centered at x,y = 0,0) as  x^2 + y^2 = radius^2,   then
1203   // (x0 + t*(x1-z0))^2 + (y0+t*(y1-y0))^2 = radius^2
1204   // (x0^2 + y0^2 - radius^2) + (2x0*(x1-x0) + 2y0*(y1-y0))*t +  ((x1-x0)^2 + (y1-y0)^2)*t^2 = 0 = C + B*t + A*t^2
1205   // quadratic with:  A = (x1-x0)^2+(y1-y0)^2 ;  B = 2x0*(x1-x0) + 2y0*(y1-y0);  C = x0^2 + y0^2 - radius^2
1206   // solution: t = (-B +/- std::sqrt(B^2 - 4*A*C)) / (2*A)
1207
1208   float A = (x[1] - x[0]) * (x[1] - x[0]) + (y[1] - y[0]) * (y[1] - y[0]);
1209   float B = 2.0 * x[0] * (x[1] - x[0]) + 2.0 * y[0] * (y[1] - y[0]);
1210   float C = x[0] * x[0] + y[0] * y[0] - radius * radius;
1211   float tup = (-B + std::sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
1212   float tdn = (-B - std::sqrt(B * B - 4.0 * A * C)) / (2.0 * A);
1213
1214   // The limits are 0 and 1, but we allow a little for floating point precision
1215   float t;
1216   if (tdn >= -0.0e-4 && tdn <= 1.0004)
1217   {
1218     t = tdn;
1219   }
1220   else if (tup >= -0.0e-4 && tup <= 1.0004)
1221   {
1222     t = tup;
1223   }
1224   else
1225   {
1226     std::cout << PHWHERE << "   **** Oops! No valid solution for tup or tdn, tdn = " << tdn << " tup = " << tup << std::endl;
1227     std::cout << "   radius " << radius << " rbegin " << std::sqrt(x[0] * x[0] + y[0] * y[0]) << " rend " << std::sqrt(x[1] * x[1] + y[1] * y[1]) << std::endl;
1228     std::cout << "   x0 " << x[0] << " x1 " << x[1] << std::endl;
1229     std::cout << "   y0 " << y[0] << " y1 " << y[1] << std::endl;
1230     std::cout << "   z0 " << z[0] << " z1 " << z[1] << std::endl;
1231     std::cout << "   A " << A << " B " << B << " C " << C << std::endl;
1232
1233     t = -1;
1234   }
1235
1236   return t;
1237 }
1238
1239 unsigned int SvtxTruthEval::getAdcValue(double gedep)
1240 {
1241   // see TPC digitizer for algorithm
1242
1243   // drift electrons per GeV of energy deposited in the TPC
1244   double Ne_dEdx = 1.56;    // keV/cm
1245   double CF4_dEdx = 7.00;   // keV/cm
1246   double Ne_NTotal = 43;    // Number/cm
1247   double CF4_NTotal = 100;  // Number/cm
1248   double Tpc_NTot = 0.5 * Ne_NTotal + 0.5 * CF4_NTotal;
1249   double Tpc_dEdx = 0.5 * Ne_dEdx + 0.5 * CF4_dEdx;
1250   double Tpc_ElectronsPerKeV = Tpc_NTot / Tpc_dEdx;
1251   double electrons_per_gev = Tpc_ElectronsPerKeV * 1e6;
1252
1253   double gem_amplification = 1400;  // GEM output electrons per drifted electron
1254   double input_electrons = gedep * electrons_per_gev * gem_amplification;
1255
1256   // convert electrons after GEM to ADC output
1257   double ChargeToPeakVolts = 20;
1258   double ADCSignalConversionGain = ChargeToPeakVolts * 1.60e-04 * 2.4;             // 20 (or 30) mV/fC * fC/electron * scaleup factor
1259   double adc_input_voltage = input_electrons * ADCSignalConversionGain;            // mV, see comments above
1260   unsigned int adc_output = (unsigned int) (adc_input_voltage * 1024.0 / 2200.0);  // input voltage x 1024 channels over 2200 mV max range
1261   adc_output = std::min<unsigned int>(adc_output, 1023);
1262
1263   return adc_output;
1264 }