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