packages/trackreco/PHTruthClustering.cc

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