Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:18:14 

0001 // this is the new trackbase version
0002 
0003 #include "PHG4MvtxHitReco.h"
0004 
0005 #include <mvtx/CylinderGeom_Mvtx.h>
0006 #include <mvtx/CylinderGeom_MvtxHelper.h>
0007 #include <mvtx/MvtxHitPruner.h>
0008 
0009 #include <trackbase/ActsGeometry.h>
0010 #include <trackbase/MvtxDefs.h>
0011 #include <trackbase/TrkrDefs.h>
0012 #include <trackbase/TrkrHit.h>  // // make iwyu happy
0013 #include <trackbase/TrkrHitSet.h>
0014 #include <trackbase/TrkrHitSetContainer.h>  // make iwyu happy
0015 #include <trackbase/TrkrHitSetContainerv1.h>
0016 #include <trackbase/TrkrHitTruthAssoc.h>  // make iwyu happy
0017 #include <trackbase/TrkrHitTruthAssocv1.h>
0018 #include <trackbase/TrkrHitv2.h>  // for TrkrHit
0019 
0020 #include <g4tracking/TrkrTruthTrack.h>
0021 #include <g4tracking/TrkrTruthTrackContainer.h>
0022 #include <g4tracking/TrkrTruthTrackContainerv1.h>
0023 #include <trackbase/ClusHitsVerbosev1.h>
0024 #include <trackbase/TrkrClusterContainer.h>
0025 #include <trackbase/TrkrClusterContainerv4.h>
0026 #include <trackbase/TrkrClusterv4.h>
0027 
0028 #include <g4detectors/PHG4CylinderGeom.h>  // for PHG4CylinderGeom
0029 #include <g4detectors/PHG4CylinderGeomContainer.h>
0030 
0031 #include <g4main/PHG4Hit.h>
0032 #include <g4main/PHG4HitContainer.h>
0033 #include <g4main/PHG4TruthInfoContainer.h>
0034 #include <g4main/PHG4Utils.h>
0035 
0036 #include <cdbobjects/CDBTTree.h>
0037 #include <ffamodules/CDBInterface.h>  // for accessing the MVTX hot pixel file from the CDB
0038 #include <ffarawobjects/MvtxRawEvtHeader.h>
0039 
0040 #include <fun4all/Fun4AllReturnCodes.h>
0041 #include <fun4all/SubsysReco.h>  // for SubsysReco
0042 
0043 #include <phool/PHCompositeNode.h>
0044 #include <phool/PHIODataNode.h>
0045 #include <phool/PHNode.h>  // for PHNode
0046 #include <phool/PHNodeIterator.h>
0047 #include <phool/PHObject.h>      // for PHObject
0048 #include <phool/PHRandomSeed.h>  // for PHRandomSeed
0049 #include <phool/getClass.h>
0050 #include <phool/phool.h>  // for PHWHERE
0051 
0052 #include <G4SystemOfUnits.hh>  // for microsecond
0053 
0054 #include <TVector3.h>  // for TVector3, ope...
0055 
0056 #include <boost/format.hpp>
0057 
0058 #include <cassert>  // for assert
0059 #include <cmath>
0060 #include <cstdlib>
0061 #include <iostream>
0062 #include <memory>  // for allocator_tra...
0063 #include <set>     // for vector
0064 #include <vector>  // for vector
0065 
0066 // New headers I added
0067 
0068 PHG4MvtxHitReco::PHG4MvtxHitReco(const std::string& name, const std::string& detector)
0069   : SubsysReco(name)
0070   , PHParameterInterface(name)
0071   , m_detector(detector)
0072   , m_tmin(-5000.)
0073   , m_tmax(5000.)
0074   , m_strobe_width(5.)
0075   , m_strobe_separation(0.)
0076   , m_truth_hits{new TrkrHitSetContainerv1}
0077 {
0078   if (Verbosity())
0079   {
0080     std::cout << "Creating PHG4MvtxHitReco for detector = " << detector << std::endl;
0081   }
0082   // initialize rng
0083   const uint seed = PHRandomSeed();
0084   m_rng.reset(gsl_rng_alloc(gsl_rng_mt19937));
0085   gsl_rng_set(m_rng.get(), seed);
0086 
0087   InitializeParameters();
0088 }
0089 
0090 int PHG4MvtxHitReco::InitRun(PHCompositeNode* topNode)
0091 {
0092   UpdateParametersWithMacro();
0093 
0094   m_tmin = get_double_param("mvtx_tmin");
0095   m_tmax = get_double_param("mvtx_tmax");
0096   m_strobe_width = get_double_param("mvtx_strobe_width");
0097   m_strobe_separation = get_double_param("mvtx_strobe_separation");
0098   m_in_sphenix_srdo = (bool) get_int_param("mvtx_in_sphenix_srdo");
0099 
0100   m_extended_readout_time = m_tmax - m_strobe_width;
0101 
0102   // printout
0103   std::cout
0104       << "PHG4MvtxHitReco::InitRun\n"
0105       << " m_tmin: " << m_tmin << "ns, m_tmax: " << m_tmax << "ns\n"
0106       << " m_strobe_width: " << m_strobe_width << "\n"
0107       << " m_strobe_separation: " << m_strobe_separation << "\n"
0108       << " m_extended_readout_time: " << m_extended_readout_time << "\n"
0109       << " m_in_sphenix_srdo: " << (m_in_sphenix_srdo ? "true" : "false") << "\n"
0110       << std::endl;
0111 
0112   //! get DST node
0113   PHNodeIterator iter(topNode);
0114   auto dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0115   assert(dstNode);
0116 
0117   //! get detector run node
0118   auto runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
0119   assert(runNode);
0120 
0121   PHNodeIterator runiter(runNode);
0122   auto runDetNode = dynamic_cast<PHCompositeNode*>(runiter.findFirst("PHCompositeNode", m_detector));
0123   if (!runDetNode)
0124   {
0125     runDetNode = new PHCompositeNode(m_detector);
0126     runNode->addNode(runDetNode);
0127   }
0128   std::string paramNodeName = "G4CELLPARAM_" + m_detector;
0129   SaveToNodeTree(runDetNode, paramNodeName);
0130 
0131   // create hitset container if needed
0132   auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0133   if (!hitsetcontainer)
0134   {
0135     PHNodeIterator dstiter(dstNode);
0136     auto trkrnode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0137     if (!trkrnode)
0138     {
0139       trkrnode = new PHCompositeNode("TRKR");
0140       dstNode->addNode(trkrnode);
0141     }
0142 
0143     hitsetcontainer = new TrkrHitSetContainerv1;
0144     auto newNode = new PHIODataNode<PHObject>(hitsetcontainer, "TRKR_HITSET", "PHObject");
0145     trkrnode->addNode(newNode);
0146   }
0147 
0148   // create hit truth association if needed
0149   auto hittruthassoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0150   if (!hittruthassoc)
0151   {
0152     PHNodeIterator dstiter(dstNode);
0153     auto trkrnode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0154     if (!trkrnode)
0155     {
0156       trkrnode = new PHCompositeNode("TRKR");
0157       dstNode->addNode(trkrnode);
0158     }
0159 
0160     hittruthassoc = new TrkrHitTruthAssocv1;
0161     auto newNode = new PHIODataNode<PHObject>(hittruthassoc, "TRKR_HITTRUTHASSOC", "PHObject");
0162     trkrnode->addNode(newNode);
0163   }
0164 
0165   // get the nodes for the truth clustering
0166   m_truthtracks = findNode::getClass<TrkrTruthTrackContainer>(topNode, "TRKR_TRUTHTRACKCONTAINER");
0167   if (!m_truthtracks)
0168   {
0169     PHNodeIterator dstiter(dstNode);
0170     m_truthtracks = new TrkrTruthTrackContainerv1();
0171     auto newNode = new PHIODataNode<PHObject>(m_truthtracks, "TRKR_TRUTHTRACKCONTAINER", "PHObject");
0172     dstNode->addNode(newNode);
0173   }
0174 
0175   m_truthclusters = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_TRUTHCLUSTERCONTAINER");
0176   if (!m_truthclusters)
0177   {
0178     m_truthclusters = new TrkrClusterContainerv4;
0179     auto newNode = new PHIODataNode<PHObject>(m_truthclusters, "TRKR_TRUTHCLUSTERCONTAINER", "PHObject");
0180     dstNode->addNode(newNode);
0181   }
0182 
0183   m_truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0184   if (!m_truthinfo)
0185   {
0186     std::cout << PHWHERE << " PHG4TruthInfoContainer node not found on node tree" << std::endl;
0187     assert(m_truthinfo);
0188   }
0189 
0190   if (record_ClusHitsVerbose)
0191   {
0192     // get the node
0193     mClusHitsVerbose = findNode::getClass<ClusHitsVerbosev1>(topNode, "Trkr_TruthClusHitsVerbose");
0194     if (!mClusHitsVerbose)
0195     {
0196       PHNodeIterator dstiter(dstNode);
0197       auto DetNode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0198       if (!DetNode)
0199       {
0200         DetNode = new PHCompositeNode("TRKR");
0201         dstNode->addNode(DetNode);
0202       }
0203       mClusHitsVerbose = new ClusHitsVerbosev1();
0204       auto newNode = new PHIODataNode<PHObject>(mClusHitsVerbose, "Trkr_TruthClusHitsVerbose", "PHObject");
0205       DetNode->addNode(newNode);
0206     }
0207   }
0208 
0209   makePixelMask(m_deadPixelMap, "MVTX_DeadPixelMap", "TotalDeadPixels");
0210   makePixelMask(m_hotPixelMap, "MVTX_HotPixelMap", "TotalHotPixels");
0211 
0212   return Fun4AllReturnCodes::EVENT_OK;
0213 }
0214 
0215 int PHG4MvtxHitReco::process_event(PHCompositeNode* topNode)
0216 {
0217   ActsGeometry* tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0218   if (!tgeometry)
0219   {
0220     std::cout << "Could not locate acts geometry" << std::endl;
0221     exit(1);
0222   }
0223 
0224   // G4Hits
0225   const std::string g4hitnodename = "G4HIT_" + m_detector;
0226   auto g4hitContainer = findNode::getClass<PHG4HitContainer>(topNode, g4hitnodename);
0227   assert(g4hitContainer);
0228 
0229   // geometry
0230   const std::string geonodename = "CYLINDERGEOM_" + m_detector;
0231   auto geoNode = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename);
0232   assert(geoNode);
0233 
0234   // Get the TrkrHitSetContainer node
0235   auto trkrHitSetContainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0236   assert(trkrHitSetContainer);
0237 
0238   // Get the TrkrHitTruthAssoc node
0239   auto hitTruthAssoc = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0240   assert(hitTruthAssoc);
0241 
0242   // Generate strobe zero relative to trigger time
0243   double strobe_zero_tm_start = generate_strobe_zero_tm_start();
0244 
0245   // assumes we want the range of accepted times to be from 0 to m_extended_readout_time
0246   std::pair<double, double> alpide_pulse = generate_alpide_pulse(0.0);  // this currently just returns fixed values
0247   double clearance = 200.0;                                             // 0.2 microsecond for luck
0248   m_tmax = m_extended_readout_time + alpide_pulse.first + clearance;
0249   m_tmin = alpide_pulse.second - clearance;
0250 
0251   // The above limits will select g4hit times of 0 up to m_extended_readout_time (only) with extensions by clearance
0252   // But we really want to select all g4hit times that will be strobed, so replace clearance with something derived from
0253   // the strobe start time in future
0254 
0255   if (Verbosity() > 0)
0256   {
0257     std::cout << " m_strobe_width " << m_strobe_width << " m_strobe_separation " << m_strobe_separation << " strobe_zero_tm_start " << strobe_zero_tm_start << " m_extended_readout_time " << m_extended_readout_time << std::endl;
0258   }
0259 
0260   // loop over all of the layers in the g4hit container
0261   auto layer_range = g4hitContainer->getLayers();
0262   for (auto layer_it = layer_range.first; layer_it != layer_range.second; ++layer_it)
0263   {
0264     // get layer
0265     const auto layer = *layer_it;
0266     assert(layer < 3);
0267 
0268     // we need the geometry object for this layer
0269     auto layergeom = dynamic_cast<CylinderGeom_Mvtx*>(geoNode->GetLayerGeom(layer));
0270     assert(layergeom);
0271 
0272     // loop over the hits in this layer
0273     const PHG4HitContainer::ConstRange g4hit_range = g4hitContainer->getHits(layer);
0274 
0275     // Get some layer parameters for later use
0276     double xpixw = layergeom->get_pixel_x();
0277     double xpixw_half = xpixw / 2.0;
0278     double zpixw = layergeom->get_pixel_z();
0279     double zpixw_half = zpixw / 2.0;
0280     int maxNX = layergeom->get_NX();
0281     int maxNZ = layergeom->get_NZ();
0282 
0283     // Now loop over all g4 hits for this layer
0284     for (auto g4hit_it = g4hit_range.first; g4hit_it != g4hit_range.second; ++g4hit_it)
0285     {
0286       // get hit
0287       auto g4hit = g4hit_it->second;
0288 
0289       truthcheck_g4hit(g4hit, topNode);
0290 
0291       // cout << "From PHG4MvtxHitReco: Call hit print method: " << std::endl;
0292       if (Verbosity() > 4)
0293       {
0294         g4hit->print();
0295       }
0296 
0297       unsigned int n_replica = 1;
0298 
0299       // Function returns ns
0300       // std::pair<double, double> alpide_pulse = generate_alpide_pulse(g4hit->get_edep());
0301 
0302       double lead_edge = g4hit->get_t(0) * ns + alpide_pulse.first;
0303       double fall_edge = g4hit->get_t(1) * ns + alpide_pulse.second;
0304 
0305       if (Verbosity() > 0)
0306       {
0307         std::cout << " MvtxHitReco: t0 " << g4hit->get_t(0) << " t1 " << g4hit->get_t(1) << " lead_edge " << lead_edge
0308                   << " fall_edge " << fall_edge << " tmin " << m_tmin << " tmax " << m_tmax << std::endl;
0309       }
0310 
0311       // check that the signal occurred witin the time window 0 to extended_readout_time, discard if not
0312       if (lead_edge > m_tmax or fall_edge < m_tmin)
0313       {
0314         continue;
0315       }
0316 
0317       double t0_strobe_frame = get_strobe_frame(lead_edge, strobe_zero_tm_start);
0318       double t1_strobe_frame = get_strobe_frame(fall_edge, strobe_zero_tm_start);
0319       n_replica += t1_strobe_frame - t0_strobe_frame;
0320 
0321       if (Verbosity() > 1)
0322       {
0323         std::cout
0324             << "MVTX is in strobed timing mode\n"
0325             << "layer " << layer << " t0(ns) " << g4hit->get_t(0) << " t1(ns) " << g4hit->get_t(1) << "\n"
0326             << "strobe_zero_tm_start(us): " << strobe_zero_tm_start / microsecond << "\n"
0327             << "strobe width(us): " << m_strobe_width / microsecond << "\n"
0328             << "strobe separation(us): " << m_strobe_separation / microsecond << "\n"
0329             << "alpide pulse start(us): " << alpide_pulse.first / microsecond << "\n"
0330             << "alpide pulse end(us): " << alpide_pulse.second / microsecond << "\n"
0331             << "tm_zero_strobe_frame: " << t0_strobe_frame << "\n"
0332             << "tm_one_strobe_frame: " << t1_strobe_frame << "\n"
0333             << "number of hit replica: " << n_replica << "\n"
0334             << std::endl;
0335       }
0336 
0337       assert(n_replica >= 1);
0338 
0339       // get_property_int(const PROPERTY prop_id) const {return INT_MIN;}
0340       int stave_number = g4hit->get_property_int(PHG4Hit::prop_stave_index);
0341       int chip_number = g4hit->get_property_int(PHG4Hit::prop_chip_index);
0342 
0343       TVector3 local_in(g4hit->get_local_x(0), g4hit->get_local_y(0), g4hit->get_local_z(0));
0344       TVector3 local_out(g4hit->get_local_x(1), g4hit->get_local_y(1), g4hit->get_local_z(1));
0345       TVector3 midpoint((local_in.X() + local_out.X()) / 2.0, (local_in.Y() + local_out.Y()) / 2.0, (local_in.Z() + local_out.Z()) / 2.0);
0346 
0347       if (Verbosity() > 2)
0348       {
0349         std::cout
0350             << " world entry point position: "
0351             << g4hit->get_x(0) << " " << g4hit->get_y(0) << " " << g4hit->get_z(0) << "\n"
0352             << " world exit  point position: "
0353             << g4hit->get_x(1) << " " << g4hit->get_y(1) << " " << g4hit->get_z(1) << "\n"
0354             << " local coords of entry point from G4 "
0355             << g4hit->get_local_x(0) << " " << g4hit->get_local_y(0) << " " << g4hit->get_local_z(0)
0356             << std::endl;
0357 
0358         TVector3 world_in(g4hit->get_x(0), g4hit->get_y(0), g4hit->get_z(0));
0359         auto hskey = MvtxDefs::genHitSetKey(layer, stave_number, chip_number, 0);
0360         auto surf = tgeometry->maps().getSiliconSurface(hskey);
0361 
0362         TVector3 local_in_check = CylinderGeom_MvtxHelper::get_local_from_world_coords(surf, tgeometry, world_in);
0363         std::cout
0364             << " local coords of entry point from geom (a check) "
0365             << local_in_check.X() << " " << local_in_check.Y() << " " << local_in_check.Z() << "\n"
0366             << " local coords of exit point from G4 "
0367             << g4hit->get_local_x(1) << " " << g4hit->get_local_y(1) << " " << g4hit->get_local_z(1)
0368             << "\n"
0369             << " local coords of exit point from geom (a check) "
0370             << local_out.X() << " " << local_out.Y() << " " << local_out.Z()
0371             << std::endl;
0372       }
0373       /*
0374       if (Verbosity() > 4)
0375       {
0376         // As a check, get the positions of the hit pixels in world coordinates from the geo object
0377               auto hskey = MvtxDefs::genHitSetKey(*layer,stave_number,chip_number,strobe);
0378               auto surf = tgeometry->maps().getSiliconSurface(hskey);
0379 
0380         TVector3 location_in = layergeom->get_world_from_local_coords(surf,tgeometry, local_in);
0381         TVector3 location_out = layergeom->get_world_from_local_coords(surf,tgeometry, local_out);
0382 
0383         std::cout
0384           << std::endl
0385           << "      PHG4MvtxHitReco:  Found world entry location from geometry for  "
0386           << " stave number " << stave_number
0387           << " half stave number " << half_stave_number
0388           << " module number " << module_number
0389           << " chip number " << chip_number
0390           << std::endl
0391           << " x = " << location_in.X()
0392           << " y = " << location_in.Y()
0393           << " z = " << location_in.Z()
0394           << " radius " << sqrt(pow(location_in.X(), 2) + pow(location_in.Y(), 2))
0395           << " angle " << atan(location_in.Y() / location_in.X())
0396           << std::endl;
0397           << "     PHG4MvtxHitReco: The world entry location from G4 was "
0398           << " x = " << g4hit->get_x(0)
0399           << " y = " << g4hit->get_y(0)
0400           << " z = " << g4hit->get_z(0)
0401           << " radius " << sqrt(pow(g4hit->get_x(0), 2) + pow(g4hit->get_y(0), 2))
0402           << " angle " << atan(g4hit->get_y(0) / g4hit->get_x(0))
0403           << std::endl;
0404           << " difference in x = " << g4hit->get_x(0) - location_in.X()
0405           << " difference in y = " << g4hit->get_y(0) - location_in.Y()
0406           << " difference in z = " << g4hit->get_z(0) - location_in.Z()
0407           << " difference in radius = "
0408           << sqrt(pow(g4hit->get_x(0), 2) + pow(g4hit->get_y(0), 2)) - sqrt(pow(location_in.X(), 2) + pow(location_in.Y(), 2))
0409           << " in angle = " << atan(g4hit->get_y(0) / g4hit->get_x(0)) - atan(location_in.Y() / location_in.X())
0410           << std::endl;
0411       }
0412   */
0413       // Get the pixel number of the entry location
0414       int pixel_number_in = layergeom->get_pixel_from_local_coords(local_in);
0415       // Get the pixel number of the exit location
0416       int pixel_number_out = layergeom->get_pixel_from_local_coords(local_out);
0417 
0418       if (pixel_number_in < 0 || pixel_number_out < 0)
0419       {
0420         std::cout
0421             << "Oops!  got negative pixel number in layer " << layergeom->get_layer()
0422             << " pixel_number_in " << pixel_number_in
0423             << " pixel_number_out " << pixel_number_out
0424             << " local_in = " << local_in.X() << " " << local_in.Y() << " " << local_in.Z()
0425             << " local_out = " << local_out.X() << " " << local_out.Y() << " " << local_out.Z()
0426             << std::endl;
0427         return Fun4AllReturnCodes::ABORTEVENT;
0428       }
0429 
0430       if (Verbosity() > 3)
0431       {
0432         std::cout
0433             << "entry pixel number " << pixel_number_in
0434             << " exit pixel number " << pixel_number_out
0435             << std::endl;
0436       }
0437 
0438       std::vector<int> vpixel;
0439       std::vector<int> vxbin;
0440       std::vector<int> vzbin;
0441       std::vector<std::pair<double, double> > venergy;
0442       // double trklen = 0.0;
0443 
0444       //===================================================
0445       // OK, now we have found which sensor the hit is in, extracted the hit
0446       // position in local sensor coordinates,  and found the pixel numbers of the
0447       // entry point and exit point
0448 
0449       //====================================================
0450       // Beginning of charge sharing implementation
0451       //    Find tracklet line inside sensor
0452       //    Divide tracklet line into n segments (vary n until answer stabilizes)
0453       //    Find centroid of each segment
0454       //    Diffuse charge at each centroid
0455       //    Apportion charge between neighboring pixels
0456       //    Add the pixel energy contributions from different track segments together
0457       //====================================================
0458 
0459       TVector3 pathvec = local_in - local_out;
0460 
0461       // See figure 7.3 of the thesis by  Lucasz Maczewski (arXiv:10053.3710) for diffusion simulations in a MAPS epitaxial layer
0462       // The diffusion widths below were inspired by those plots, corresponding to where the probability drops off to 1/3 of the peak value
0463       // However note that we make the simplifying assumption that the probability distribution is flat within this diffusion width,
0464       // while in the simulation it is not
0465       // double diffusion_width_max = 35.0e-04;   // maximum diffusion radius 35 microns, in cm
0466       // double diffusion_width_min = 12.0e-04;   // minimum diffusion radius 12 microns, in cm
0467       double diffusion_width_max = 25.0e-04;  // maximum diffusion radius 35 microns, in cm
0468       double diffusion_width_min = 8.0e-04;   // minimum diffusion radius 12 microns, in cm
0469 
0470       double ydrift_max = pathvec.Y();
0471       int nsegments = 4;
0472 
0473       // we want to make a list of all pixels possibly affected by this hit
0474       // we take the entry and exit locations in local coordinates, and build
0475       // a rectangular array of pixels that encompasses both, with "nadd" pixels added all around
0476 
0477       int xbin_in = layergeom->get_pixel_X_from_pixel_number(pixel_number_in);
0478       int zbin_in = layergeom->get_pixel_Z_from_pixel_number(pixel_number_in);
0479       int xbin_out = layergeom->get_pixel_X_from_pixel_number(pixel_number_out);
0480       int zbin_out = layergeom->get_pixel_Z_from_pixel_number(pixel_number_out);
0481 
0482       int xbin_max, xbin_min;
0483       int nadd = 2;
0484       if (xbin_in > xbin_out)
0485       {
0486         xbin_max = xbin_in + nadd;
0487         xbin_min = xbin_out - nadd;
0488       }
0489       else
0490       {
0491         xbin_max = xbin_out + nadd;
0492         xbin_min = xbin_in - nadd;
0493       }
0494 
0495       int zbin_max, zbin_min;
0496       if (zbin_in > zbin_out)
0497       {
0498         zbin_max = zbin_in + nadd;
0499         zbin_min = zbin_out - nadd;
0500       }
0501       else
0502       {
0503         zbin_max = zbin_out + nadd;
0504         zbin_min = zbin_in - nadd;
0505       }
0506 
0507       // need to check that values of xbin and zbin are within the valid range
0508       // YCM: Fix pixel range: Xbin (row) 0 to 511, Zbin (col) 0 to 1023
0509       if (xbin_min < 0)
0510       {
0511         xbin_min = 0;
0512       }
0513       if (zbin_min < 0)
0514       {
0515         zbin_min = 0;
0516       }
0517       if (xbin_max >= maxNX)
0518       {
0519         xbin_max = maxNX - 1;
0520       }
0521       if (zbin_max >= maxNZ)
0522       {
0523         xbin_max = maxNZ - 1;
0524       }
0525 
0526       if (Verbosity() > 1)
0527       {
0528         std::cout << " xbin_in " << xbin_in << " xbin_out " << xbin_out << " xbin_min " << xbin_min << " xbin_max " << xbin_max << std::endl;
0529         std::cout << " zbin_in " << zbin_in << " zbin_out " << zbin_out << " zbin_min " << zbin_min << " zbin_max " << zbin_max << std::endl;
0530       }
0531 
0532       // skip this hit if it involves an unreasonable  number of pixels
0533       // this skips it if either the xbin or ybin range traversed is greater than 8 (for 8 adding two pixels at each end makes the range 12)
0534       if (xbin_max - xbin_min > 12 || zbin_max - zbin_min > 12)
0535       {
0536         continue;
0537       }
0538 
0539       // this hit is skipped earlier if this dimensioning would be exceeded
0540       double pixenergy[12][12] = {};  // init to 0
0541       double pixeion[12][12] = {};    // init to 0
0542 
0543       // Loop over track segments and diffuse charge at each segment location, collect energy in pixels
0544       for (int i = 0; i < nsegments; i++)
0545       {
0546         // Find the tracklet segment location
0547         // If there are n segments of equal length, we want 2*n intervals
0548         // The 1st segment is centered at interval 1, the 2nd at interval 3, the nth at interval 2n -1
0549         double interval = 2 * (double) i + 1;
0550         double frac = interval / (double) (2 * nsegments);
0551         TVector3 segvec(pathvec.X() * frac, pathvec.Y() * frac, pathvec.Z() * frac);
0552         segvec = segvec + local_out;
0553 
0554         //  Find the distance to the back of the sensor from the segment location
0555         // That projection changes only the value of y
0556         double ydrift = segvec.Y() - local_out.Y();
0557 
0558         // Caculate the charge diffusion over this drift distance
0559         // increases from diffusion width_min to diffusion_width_max
0560         double ydiffusion_radius = diffusion_width_min + (ydrift / ydrift_max) * (diffusion_width_max - diffusion_width_min);
0561 
0562         if (Verbosity() > 5)
0563         {
0564           std::cout
0565               << " segment " << i
0566               << " interval " << interval
0567               << " frac " << frac
0568               << " local_in.X " << local_in.X()
0569               << " local_in.Z " << local_in.Z()
0570               << " local_in.Y " << local_in.Y()
0571               << " pathvec.X " << pathvec.X()
0572               << " pathvec.Z " << pathvec.Z()
0573               << " pathvec.Y " << pathvec.Y()
0574               << " segvec.X " << segvec.X()
0575               << " segvec.Z " << segvec.Z()
0576               << " segvec.Y " << segvec.Y()
0577               << " ydrift " << ydrift
0578               << " ydrift_max " << ydrift_max
0579               << " ydiffusion_radius " << ydiffusion_radius
0580               << std::endl;
0581         }
0582         // Now find the area of overlap of the diffusion circle with each pixel and apportion the energy
0583         for (int ix = xbin_min; ix <= xbin_max; ix++)
0584         {
0585           for (int iz = zbin_min; iz <= zbin_max; iz++)
0586           {
0587             // Find the pixel corners for this pixel number
0588             int pixnum = layergeom->get_pixel_number_from_xbin_zbin(ix, iz);
0589 
0590             if (pixnum < 0)
0591             {
0592               std::cout
0593                   << " pixnum < 0 , pixnum = " << pixnum << "\n"
0594                   << " ix " << ix << " iz " << iz << "\n"
0595                   << " xbin_min " << xbin_min << " zbin_min " << zbin_min << "\n"
0596                   << " xbin_max " << xbin_max << " zbin_max " << zbin_max << "\n"
0597                   << " maxNX " << maxNX << " maxNZ " << maxNZ
0598                   << std::endl;
0599             }
0600 
0601             TVector3 tmp = layergeom->get_local_coords_from_pixel(pixnum);
0602             // note that (x1,z1) is the top left corner, (x2,z2) is the bottom right corner of the pixel - circle_rectangle_intersection expects this ordering
0603             double x1 = tmp.X() - xpixw_half;
0604             double z1 = tmp.Z() + zpixw_half;
0605             double x2 = tmp.X() + xpixw_half;
0606             double z2 = tmp.Z() - zpixw_half;
0607 
0608             // here segvec.X and segvec.Z are the center of the circle, and diffusion_radius is the circle radius
0609             // circle_rectangle_intersection returns the overlap area of the circle and the pixel. It is very fast if there is no overlap.
0610             double pixarea_frac = PHG4Utils::circle_rectangle_intersection(x1, z1, x2, z2, segvec.X(), segvec.Z(), ydiffusion_radius) / (M_PI * pow(ydiffusion_radius, 2));
0611             // assume that the energy is deposited uniformly along the tracklet length, so that this segment gets the fraction 1/nsegments of the energy
0612             pixenergy[ix - xbin_min][iz - zbin_min] += pixarea_frac * g4hit->get_edep() / (float) nsegments;
0613             if (g4hit->has_property(PHG4Hit::prop_eion))
0614             {
0615               pixeion[ix - xbin_min][iz - zbin_min] += pixarea_frac * g4hit->get_eion() / (float) nsegments;
0616             }
0617             if (Verbosity() > 5)
0618             {
0619               std::cout
0620                   << "    pixnum " << pixnum << " xbin " << ix << " zbin " << iz
0621                   << " pixel_area fraction of circle " << pixarea_frac << " accumulated pixel energy " << pixenergy[ix - xbin_min][iz - zbin_min]
0622                   << std::endl;
0623             }
0624           }
0625         }
0626       }  // end loop over segments
0627 
0628       // now we have the energy deposited in each pixel, summed over all tracklet segments. We make a vector of all pixels with non-zero energy deposited
0629       for (int ix = xbin_min; ix <= xbin_max; ix++)
0630       {
0631         for (int iz = zbin_min; iz <= zbin_max; iz++)
0632         {
0633           if (pixenergy[ix - xbin_min][iz - zbin_min] > 0.0)
0634           {
0635             int pixnum = layergeom->get_pixel_number_from_xbin_zbin(ix, iz);
0636             vpixel.push_back(pixnum);
0637             vxbin.push_back(ix);
0638             vzbin.push_back(iz);
0639             std::pair<double, double> tmppair = std::make_pair(pixenergy[ix - xbin_min][iz - zbin_min], pixeion[ix - xbin_min][iz - zbin_min]);
0640             venergy.push_back(tmppair);
0641             if (Verbosity() > 1)
0642             {
0643               std::cout
0644                   << " Added pixel number " << pixnum << " xbin " << ix
0645                   << " zbin " << iz << " to vectors with energy " << pixenergy[ix - xbin_min][iz - zbin_min]
0646                   << std::endl;
0647             }
0648           }
0649         }
0650       }
0651 
0652       //===================================
0653       // End of charge sharing implementation
0654       //===================================
0655 
0656       // loop over all fired cells for this g4hit and add them to the TrkrHitSet
0657 
0658       for (unsigned int i1 = 0; i1 < vpixel.size(); i1++)  // loop over all fired cells
0659       {
0660         // This is the new storage object version
0661         //====================================
0662         for (unsigned int i_rep = 0; i_rep < n_replica; i_rep++)
0663         {
0664           int strobe = t0_strobe_frame + i_rep;
0665           // to fit in a 5 bit field in the hitsetkey [-16,15]
0666           if (strobe < -16)
0667           {
0668             strobe = -16;
0669           }
0670           if (strobe >= 16)
0671           {
0672             strobe = 15;
0673           }
0674 
0675           // We need to create the TrkrHitSet if not already made - each TrkrHitSet should correspond to a chip for the Mvtx
0676           TrkrDefs::hitsetkey hitsetkey = MvtxDefs::genHitSetKey(layer, stave_number, chip_number, strobe);
0677           // Use existing hitset or add new one if needed
0678           TrkrHitSetContainer::Iterator hitsetit = trkrHitSetContainer->findOrAddHitSet(hitsetkey);
0679 
0680           // generate the key for this hit
0681           TrkrDefs::hitkey hitkey = MvtxDefs::genHitKey(vzbin[i1], vxbin[i1]);
0682           // See if this hit already exists
0683           TrkrHit* hit = nullptr;
0684           hit = hitsetit->second->getHit(hitkey);
0685 
0686           if (hit)
0687           {
0688             if (Verbosity() > 0)
0689             {
0690               std::cout << PHWHERE << "::" << __func__
0691                         << " - duplicated hit, hitsetkey: " << hitsetkey
0692                         << " hitkey: " << hitkey << std::endl;
0693             }
0694             continue;
0695           }
0696           const TrkrDefs::hitsetkey hitsetkeymask = MvtxDefs::genHitSetKey(layer, stave_number, chip_number, 0);
0697 
0698           // Regardless of whether the hit should be masked, add the energy to the truth hit
0699           double hitenergy = venergy[i1].first * TrkrDefs::MvtxEnergyScaleup;
0700           addtruthhitset(hitsetkey, hitkey, hitenergy);
0701 
0702           if ((std::find(m_deadPixelMap.begin(), m_deadPixelMap.end(), std::make_pair(hitsetkeymask, hitkey)) == m_deadPixelMap.end()) && (std::find(m_hotPixelMap.begin(), m_hotPixelMap.end(), std::make_pair(hitsetkeymask, hitkey)) == m_hotPixelMap.end()))
0703           {
0704             // create hit and insert in hitset
0705             hit = new TrkrHitv2();
0706 
0707             hit->addEnergy(hitenergy);
0708             hitsetit->second->addHitSpecificKey(hitkey, hit);
0709           }
0710           else
0711           {
0712             continue;
0713           }
0714 
0715           addtruthhitset(hitsetkey, hitkey, hitenergy);
0716 
0717           if (Verbosity() > 0)
0718           {
0719             std::cout << "Layer: " << layer << ", Stave: " << (uint16_t) MvtxDefs::getStaveId(hitsetkey) << ", Chip: " << (uint16_t) MvtxDefs::getChipId(hitsetkey) << ", Row: " << (uint16_t) MvtxDefs::getRow(hitkey) << ", Col: " << (uint16_t) MvtxDefs::getCol(hitkey) << ", Strobe: " << (int) MvtxDefs::getStrobeId(hitsetkey) << ", added hit " << hitkey << " to hitset " << hitsetkey << " with energy " << hit->getEnergy() / TrkrDefs::MvtxEnergyScaleup << std::endl;
0720           }
0721 
0722           // now we update the TrkrHitTruthAssoc map - the map contains <hitsetkey, std::pair <hitkey, g4hitkey> >
0723           // There is only one TrkrHit per pixel, but there may be multiple g4hits
0724           // How do we know how much energy from PHG4Hit went into TrkrHit? We don't, have to sort it out in evaluator to save memory
0725 
0726           // we set the strobe ID to zero in the hitsetkey
0727           // we use the findOrAdd method to keep from adding identical entries
0728           TrkrDefs::hitsetkey bare_hitsetkey = zero_strobe_bits(hitsetkey);
0729           hitTruthAssoc->findOrAddAssoc(bare_hitsetkey, hitkey, g4hit_it->first);
0730         }
0731       }  // end loop over hit cells
0732     }    // end loop over g4hits for this layer
0733 
0734   }  // end loop over layers
0735 
0736   // print the list of entries in the association table
0737   if (Verbosity() > 0)
0738   {
0739     std::cout << "From PHG4MvtxHitReco: " << std::endl;
0740     hitTruthAssoc->identify();
0741   }
0742 
0743   // spit out the truth clusters
0744 
0745   end_event_truthcluster(topNode);
0746 
0747   if (Verbosity() > 3)
0748   {
0749     int nclusprint = -1;
0750     std::cout << PHWHERE << ": content of clusters " << std::endl;
0751     auto& tmap = m_truthtracks->getMap();
0752     std::cout << " Number of tracks: " << tmap.size() << std::endl;
0753     for (auto& _pair : tmap)
0754     {
0755       auto& track = _pair.second;
0756       std::cout << (boost::format("id(%2i) phi:eta:pt(%5.2f:%5.2f:%5.2f) nclusters(") % (int) track->getTrackid() % track->getPhi() % track->getPseudoRapidity() % track->getPt()).str() << track->getClusters().size() << ") " << std::endl;
0757       int nclus = 0;
0758       for (auto cluskey : track->getClusters())
0759       {
0760         std::cout << " "
0761                   << ((int) TrkrDefs::getHitSetKeyFromClusKey(cluskey)) << ":index(" << ((int) TrkrDefs::getClusIndex(cluskey)) << ")" << std::endl;
0762         ++nclus;
0763         if (nclusprint > 0 && nclus >= nclusprint)
0764         {
0765           std::cout << " ... ";
0766           break;
0767         }
0768       }
0769     }
0770     std::cout << PHWHERE << " ----- end of clusters " << std::endl;
0771   }
0772 
0773   if (m_is_emb)
0774   {
0775     cluster_truthhits(topNode);  // the last track was truth -- cluster it
0776     prior_g4hit = nullptr;
0777   }
0778 
0779   return Fun4AllReturnCodes::EVENT_OK;
0780 }
0781 
0782 std::pair<double, double> PHG4MvtxHitReco::generate_alpide_pulse(const double energy_deposited)
0783 {
0784   // We need to translate energy deposited to num/ electrons released
0785   if (Verbosity() > 2)
0786   {
0787     std::cout << "energy_deposited: " << energy_deposited << std::endl;
0788   }
0789   // int silicon_band_gap = 1.12; //Band gap energy in eV
0790   // int Q_in = rand() % 5000 + 50;
0791   // int clipping_point = 110;
0792   // double ToT_start = Q_in < 200 ? 395.85*exp(-0.5*pow((Q_in+851.43)/286.91, 2)) : 0.5;
0793   // double ToT_end = Q_in < clipping_point ? 5.90*exp(-0.5*pow((Q_in-99.86)/54.80, 2)) : 5.8 - 6.4e-4 * Q_in;
0794 
0795   // return make_pair(ToT_start*1e3, ToT_end*1e3);
0796   //  Using constant alpide pulse length
0797   return std::make_pair<double, double>(1.5 * microsecond, 5.9 * microsecond);
0798 }
0799 
0800 double PHG4MvtxHitReco::generate_strobe_zero_tm_start()
0801 {
0802   return -1. * gsl_rng_uniform_pos(m_rng.get()) * (m_strobe_separation + m_strobe_width);
0803 }
0804 
0805 int PHG4MvtxHitReco::get_strobe_frame(double alpide_time, double strobe_zero_tm_start)
0806 {
0807   int strobe_frame = int((alpide_time - strobe_zero_tm_start) / (m_strobe_width + m_strobe_separation));
0808   strobe_frame += (alpide_time < strobe_zero_tm_start) ? -1 : 0;
0809   return strobe_frame;
0810 }
0811 
0812 void PHG4MvtxHitReco::set_timing_window(const int detid, const double tmin, const double tmax)
0813 {
0814   if (false)
0815   {
0816     std::cout
0817         << "PHG4MvtxHitReco: Set Mvtx timing window parameters from macro for layer = "
0818         << detid << " to tmin = " << tmin << " tmax = " << tmax
0819         << std::endl;
0820   }
0821 
0822   return;
0823 }
0824 
0825 void PHG4MvtxHitReco::SetDefaultParameters()
0826 {
0827   // cout << "PHG4MvtxHitReco: Setting Mvtx timing window defaults to tmin = -5000 and  tmax = 5000 ns" << std::endl;
0828   set_default_double_param("mvtx_tmin", -5000);
0829   set_default_double_param("mvtx_tmax", 5000);
0830   set_default_double_param("mvtx_strobe_width", 5 * microsecond);
0831   set_default_double_param("mvtx_strobe_separation", 0.);
0832   set_default_int_param("mvtx_in_sphenix_srdo", (int) false);
0833   return;
0834 }
0835 
0836 TrkrDefs::hitsetkey PHG4MvtxHitReco::zero_strobe_bits(TrkrDefs::hitsetkey hitsetkey)
0837 {
0838   unsigned int layer = TrkrDefs::getLayer(hitsetkey);
0839   unsigned int stave = MvtxDefs::getStaveId(hitsetkey);
0840   unsigned int chip = MvtxDefs::getChipId(hitsetkey);
0841   TrkrDefs::hitsetkey bare_hitsetkey = MvtxDefs::genHitSetKey(layer, stave, chip, 0);
0842 
0843   return bare_hitsetkey;
0844 }
0845 
0846 PHG4MvtxHitReco::~PHG4MvtxHitReco()
0847 {
0848   delete m_truth_hits;
0849 }
0850 
0851 void PHG4MvtxHitReco::truthcheck_g4hit(PHG4Hit* g4hit, PHCompositeNode* topNode)
0852 {
0853   int new_trkid = (g4hit == nullptr) ? -1 : g4hit->get_trkid();
0854   bool is_new_track = (new_trkid != m_trkid);
0855   if (Verbosity() > 5)
0856   {
0857     std::cout << PHWHERE << std::endl
0858               << " -> Checking status of PHG4Hit. Track id(" << new_trkid << ")" << std::endl;
0859   }
0860   if (!is_new_track)
0861   {
0862     /* FIXME */
0863     /* std::cout << " FIXME PEAR Checking g4hit : " << g4hit->get_x(0) << " " */
0864     /*   << g4hit->get_y(0) << " " << g4hit->get_z(0) */
0865     /*   << "  trackid("<<g4hit->get_trkid() << ") " << std::endl; */
0866     if (m_is_emb)
0867     {
0868       /* std::cout << " FIXME Checking MVTX g4hit : " << g4hit->get_x(0) << " " */
0869       /*   << g4hit->get_y(0) << " " << g4hit->get_z(0) */
0870       /*   << "  trackid("<<g4hit->get_trkid() << ") " << std::endl; */
0871       if (prior_g4hit != nullptr && (std::abs(prior_g4hit->get_x(0) - g4hit->get_x(0)) > max_g4hitstep || std::abs(prior_g4hit->get_y(0) - g4hit->get_y(0)) > max_g4hitstep))
0872       {
0873         // this is a looper track -- cluster hits up to this point already
0874         cluster_truthhits(topNode);
0875       }
0876       prior_g4hit = g4hit;
0877     }
0878     return;
0879   }
0880   // <- STATUS: this is a new track
0881   if (Verbosity() > 2)
0882   {
0883     std::cout << PHWHERE << std::endl
0884               << " -> Found new embedded track with id: " << new_trkid << std::endl;
0885   }
0886   if (m_is_emb)
0887   {
0888     cluster_truthhits(topNode);  // cluster m_truth_hits and add m_current_track
0889     m_current_track = nullptr;
0890     prior_g4hit = nullptr;
0891   }
0892   m_trkid = new_trkid;
0893   m_is_emb = m_truthinfo->isEmbeded(m_trkid);
0894   if (m_is_emb)
0895   {
0896     m_current_track = m_truthtracks->getTruthTrack(m_trkid, m_truthinfo);
0897     prior_g4hit = g4hit;
0898   }
0899 }
0900 
0901 void PHG4MvtxHitReco::end_event_truthcluster(PHCompositeNode* topNode)
0902 {
0903   if (m_is_emb)
0904   {
0905     cluster_truthhits(topNode);  // cluster m_truth_hits and add m_current_track
0906     m_current_track = nullptr;
0907     m_trkid = -1;
0908     m_is_emb = false;
0909   }
0910   m_hitsetkey_cnt.clear();
0911 }
0912 
0913 void PHG4MvtxHitReco::addtruthhitset(
0914     TrkrDefs::hitsetkey hitsetkey,
0915     TrkrDefs::hitkey hitkey,
0916     float neffelectrons)
0917 {
0918   if (!m_is_emb)
0919   {
0920     return;
0921   }
0922   TrkrHitSetContainer::Iterator hitsetit = m_truth_hits->findOrAddHitSet(hitsetkey);
0923   // See if this hit already exists
0924   TrkrHit* hit = nullptr;
0925   hit = hitsetit->second->getHit(hitkey);
0926   if (!hit)
0927   {
0928     // create a new one
0929     hit = new TrkrHitv2();
0930     hitsetit->second->addHitSpecificKey(hitkey, hit);
0931   }
0932   // Either way, add the energy to it  -- adc values will be added at digitization
0933   hit->addEnergy(neffelectrons);
0934 }
0935 
0936 void PHG4MvtxHitReco::cluster_truthhits(PHCompositeNode* topNode)
0937 {
0938   // clusterize the mvtx hits in m_truth_hits, put them in m_truthclusters, and put the id's in m_current_track
0939   // ----------------------------------------------------------------------------------------------------
0940   // Digitization -- simplified from g4mvtx/PHG4MvtxDigitizer --
0941   // ----------------------------------------------------------------------------------------------------
0942 
0943   /* // We want all hitsets for the Mvtx */
0944   /* TrkrHitSetContainer::ConstRange hitset_range = m_truth_hits->getHitSets(TrkrDefs::TrkrId::mvtxId); */
0945   /* for (TrkrHitSetContainer::ConstIterator hitset_iter = hitset_range.first; */
0946   /*      hitset_iter != hitset_range.second; */
0947   /*      ++hitset_iter) */
0948   /* { */
0949   /*   // we have an iterator to one TrkrHitSet for the mvtx from the trkrHitSetContainer */
0950   /*   // get the hitset key so we can find the layer */
0951   /*   TrkrDefs::hitsetkey hitsetkey = hitset_iter->first; */
0952   /*   int layer = TrkrDefs::getLayer(hitsetkey); */
0953   /*   // FIXME */
0954   /*   /1* if (Verbosity() > 1) std::cout << "PHG4MvtxDigitizer: found hitset with key: " << hitsetkey << " in layer " << layer << std::endl; *1/ */
0955   /*   std::cout << "PHG4MvtxDigitizer: found hitset with key: " << hitsetkey << " in layer " << layer << std::endl; */
0956 
0957   /*   // get all of the hits from this hitset */
0958   /*   TrkrHitSet *hitset = hitset_iter->second; */
0959   /*   TrkrHitSet::ConstRange hit_range = hitset->getHits(); */
0960   /*   std::set<TrkrDefs::hitkey> hits_rm; */
0961   /*   for (TrkrHitSet::ConstIterator hit_iter = hit_range.first; */
0962   /*        hit_iter != hit_range.second; */
0963   /*        ++hit_iter) */
0964   /*   { */
0965   /*     TrkrHit *hit = hit_iter->second; */
0966 
0967   /*     // Convert the signal value to an ADC value and write that to the hit */
0968   /*     // Unsigned int adc = hit->getEnergy() / (TrkrDefs::MvtxEnergyScaleup *_energy_scale[layer]); */
0969   /*     if (Verbosity() > 0) */
0970   /*       std::cout << "    PHG4MvtxDigitizer: found hit with key: " << hit_iter->first << " and signal " << hit->getEnergy() / TrkrDefs::MvtxEnergyScaleup << " in layer " << layer << std::endl; */
0971   /*     // Remove the hits with energy under threshold */
0972   /*     bool rm_hit = false; */
0973   /*     // FIXME: */
0974   /*     double _energy_threshold = 0.; // FIXME */
0975   /*     const double dummy_pixel_thickness = 0.001; */
0976   /*     const double mip_e = 0.003876; */
0977   /*     double _energy_scale =  mip_e * dummy_pixel_thickness; // FIXME: note: this doesn't actually */
0978   /*                                                            // matter either here or for the Svtx Tracks -- the energy is digital -- either the hit is there or it isn't */
0979   /*     double _max_adc = 255; */
0980   /*     if ((hit->getEnergy() / TrkrDefs::MvtxEnergyScaleup) < _energy_threshold) */
0981   /*     { */
0982   /*       if (Verbosity() > 0) std::cout << "         remove hit, below energy threshold of " << _energy_threshold << std::endl; */
0983   /*       rm_hit = true; */
0984   /*     } */
0985   /*     unsigned short adc = (unsigned short) (hit->getEnergy() / (TrkrDefs::MvtxEnergyScaleup * _energy_scale)); */
0986   /*     if (adc > _max_adc) adc = _max_adc; */
0987   /*     hit->setAdc(adc); */
0988 
0989   /* if (rm_hit) hits_rm.insert(hit_iter->first); */
0990   /*   } */
0991 
0992   /*   for (const auto &key : hits_rm) */
0993   /*   { */
0994   /*     if (Verbosity() > 0) std::cout << "    PHG4MvtxDigitizer: remove hit with key: " << key << std::endl; */
0995   /*     hitset->removeHit(key); */
0996   /*   } */
0997   /* } */
0998 
0999   // ----------------------------------------------------------------------------------------------------
1000   // Prune hits -- simplified from mvtx/MvtxHitReco
1001   // ----------------------------------------------------------------------------------------------------
1002   std::multimap<TrkrDefs::hitsetkey, TrkrDefs::hitsetkey> hitset_multimap;  // will map (bare hitset, hitset with strobe)
1003   std::set<TrkrDefs::hitsetkey> bare_hitset_set;                            // list of all physical sensor hitsetkeys (i.e. with strobe set to zero)
1004 
1005   TrkrHitSetContainer::ConstRange hitsetrange = m_truth_hits->getHitSets(TrkrDefs::TrkrId::mvtxId);  // actually m_truth_hits can only have MVTX hits at this point...
1006   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
1007        hitsetitr != hitsetrange.second;
1008        ++hitsetitr)
1009   {
1010     auto hitsetkey = hitsetitr->first;
1011 
1012     // get the hitsetkey value for strobe 0
1013     unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
1014     unsigned int stave = MvtxDefs::getStaveId(hitsetitr->first);
1015     unsigned int chip = MvtxDefs::getChipId(hitsetitr->first);
1016     auto bare_hitsetkey = MvtxDefs::genHitSetKey(layer, stave, chip, 0);
1017 
1018     hitset_multimap.insert(std::make_pair(bare_hitsetkey, hitsetkey));
1019     bare_hitset_set.insert(bare_hitsetkey);
1020 
1021     if (Verbosity() > 0)
1022     {
1023       std::cout << " found hitsetkey " << hitsetkey << " for bare_hitsetkey " << bare_hitsetkey << std::endl;
1024     }
1025   }
1026 
1027   // Now consolidate all hits into the hitset with strobe 0, and delete the other hitsets
1028   //==============================================================
1029   for (unsigned int bare_hitsetkey : bare_hitset_set)
1030   {
1031     TrkrHitSet* bare_hitset = (m_truth_hits->findOrAddHitSet(bare_hitsetkey))->second;
1032     if (Verbosity() > 0)
1033     {
1034       std::cout << "         bare_hitset " << bare_hitsetkey << " initially has " << bare_hitset->size() << " hits " << std::endl;
1035     }
1036 
1037     auto bare_hitsetrange = hitset_multimap.equal_range(bare_hitsetkey);
1038     for (auto it = bare_hitsetrange.first; it != bare_hitsetrange.second; ++it)
1039     {
1040       auto hitsetkey = it->second;
1041 
1042       int strobe = MvtxDefs::getStrobeId(hitsetkey);
1043       if (strobe != 0)
1044       {
1045         if (Verbosity() > 0)
1046         {
1047           std::cout << "            process hitsetkey " << hitsetkey << " for bare_hitsetkey " << bare_hitsetkey << std::endl;
1048         }
1049 
1050         // copy all hits to the hitset with strobe 0
1051         TrkrHitSet* hitset = m_truth_hits->findHitSet(hitsetkey);
1052 
1053         if (Verbosity() > 0)
1054         {
1055           std::cout << "                hitsetkey " << hitsetkey << " has strobe " << strobe << " and has " << hitset->size() << " hits,  so copy it" << std::endl;
1056         }
1057 
1058         TrkrHitSet::ConstRange hitrangei = hitset->getHits();
1059         for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
1060              hitr != hitrangei.second;
1061              ++hitr)
1062         {
1063           auto hitkey = hitr->first;
1064           if (Verbosity() > 0)
1065           {
1066             std::cout << "                 found hitkey " << hitkey << std::endl;
1067           }
1068           // if it is already there, leave it alone, this is a duplicate hit
1069           auto tmp_hit = bare_hitset->getHit(hitkey);
1070           if (tmp_hit)
1071           {
1072             if (Verbosity() > 0)
1073             {
1074               std::cout << "                          hitkey " << hitkey << " is already in bare hitsest, do not copy" << std::endl;
1075             }
1076             continue;
1077           }
1078 
1079           // otherwise copy the hit over
1080           if (Verbosity() > 0)
1081           {
1082             std::cout << "                          copying over hitkey " << hitkey << std::endl;
1083           }
1084           auto old_hit = hitr->second;
1085           TrkrHit* new_hit = new TrkrHitv2();
1086           new_hit->setAdc(old_hit->getAdc());
1087           bare_hitset->addHitSpecificKey(hitkey, new_hit);
1088         }
1089 
1090         // all hits are copied over to the strobe zero hitset, remove this hitset
1091         m_truth_hits->removeHitSet(hitsetkey);
1092       }
1093     }
1094   }
1095 
1096   // ----------------------------------------------------------------------------------------------------
1097   // cluster hits -- simplified from mvtx/MvtxClusterizer
1098   // ----------------------------------------------------------------------------------------------------
1099   PHG4CylinderGeomContainer* geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
1100   if (!geom_container)
1101   {
1102     std::cout << PHWHERE << std::endl;
1103     std::cout << "WARNING: cannot find the geometry CYLINDERGEOM_MVTX" << std::endl;
1104     m_truth_hits->Reset();
1105     return;
1106   }
1107 
1108   //-----------
1109   // Clustering
1110   //-----------
1111 
1112   // loop over each MvtxHitSet object (chip)
1113   hitsetrange = m_truth_hits->getHitSets(TrkrDefs::TrkrId::mvtxId);
1114   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
1115        hitsetitr != hitsetrange.second;
1116        ++hitsetitr)
1117   {                                          // hitsetitr    : pair(TrkrDefs::hitsetkey, TrkrHitSet>;   TrkrHitSet : map <HitKey, TrkrHit>
1118     TrkrHitSet* hitset = hitsetitr->second;  // hitset : map <TrkrDefs::hitkey, TrkrHit>
1119 
1120     /* if (true) */
1121     if (Verbosity() > 0)
1122     {
1123       unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
1124       unsigned int stave = MvtxDefs::getStaveId(hitsetitr->first);
1125       unsigned int chip = MvtxDefs::getChipId(hitsetitr->first);
1126       unsigned int strobe = MvtxDefs::getStrobeId(hitsetitr->first);
1127       std::cout << "MvtxClusterizer found hitsetkey " << hitsetitr->first << " layer " << layer << " stave " << stave << " chip " << chip << " strobe " << strobe << std::endl;
1128     }
1129 
1130     if (Verbosity() > 2)
1131     {
1132       hitset->identify();
1133     }
1134 
1135     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
1136 
1137     auto hitsetkey = hitset->getHitSetKey();
1138     auto ckey = TrkrDefs::genClusKey(hitsetkey, 0);  // there is only one cluster made per cluskey
1139 
1140     // determine the size of the cluster in phi and z
1141     std::set<int> phibins;
1142     std::set<int> zbins;
1143 
1144     // determine the cluster position...
1145     double locxsum = 0.;
1146     double loczsum = 0.;
1147 
1148     double locclusx = NAN;
1149     double locclusz = NAN;
1150 
1151     // we need the geometry object for this layer to get the global positions
1152     int layer = TrkrDefs::getLayer(ckey);
1153     auto layergeom = dynamic_cast<CylinderGeom_Mvtx*>(geom_container->GetLayerGeom(layer));
1154     if (!layergeom)
1155     {
1156       exit(1);
1157     }
1158 
1159     // make a tunable threshold for energy in a given hit
1160     //  -- percentage of someting? (total cluster energy)
1161     double sum_adc{0.};
1162     for (auto ihit = hitrangei.first; ihit != hitrangei.second; ++ihit)
1163     {
1164       sum_adc += ihit->second->getAdc();
1165     }
1166     const double threshold = sum_adc * m_pixel_thresholdrat;       // FIXME -- tune this as needed
1167     std::map<int, unsigned int> m_iphi, m_it, m_iphiCut, m_itCut;  // FIXME
1168                                                                    /* const unsigned int npixels = std::distance( hitrangei.first, hitrangei.second ); */
1169     // to tune this parameter: run a bunch of tracks and compare truth sizes and reco sizes,
1170     // should come out the same
1171     double npixels{0.};
1172     for (auto ihit = hitrangei.first; ihit != hitrangei.second; ++ihit)
1173     {
1174       const auto adc = ihit->second->getAdc();
1175       int col = MvtxDefs::getCol(ihit->first);
1176       int row = MvtxDefs::getRow(ihit->first);
1177 
1178       if (mClusHitsVerbose)
1179       {
1180         std::map<int, unsigned int>& m_phi = (adc < threshold) ? m_iphiCut : m_iphi;
1181         std::map<int, unsigned int>& m_z = (adc < threshold) ? m_itCut : m_it;
1182 
1183         auto pnew = m_phi.try_emplace(row, adc);
1184         if (!pnew.second)
1185         {
1186           pnew.first->second += adc;
1187         }
1188 
1189         pnew = m_z.try_emplace(col, adc);
1190         if (!pnew.second)
1191         {
1192           pnew.first->second += adc;
1193         }
1194       }
1195       if (adc < threshold)
1196       {
1197         continue;
1198       }
1199 
1200       // size
1201       npixels += 1.;
1202       zbins.insert(col);
1203       phibins.insert(row);
1204 
1205       // get local coordinates, in stave reference frame, for hit
1206       auto local_coords = layergeom->get_local_coords_from_pixel(row, col);
1207 
1208       /*
1209          manually offset position along y (thickness of the sensor),
1210          to account for effective hit position in the sensor, resulting from diffusion.
1211          Effective position corresponds to 1um above the middle of the sensor
1212          */
1213       local_coords.SetY(1e-4);
1214 
1215       // update cluster position
1216       locxsum += local_coords.X();
1217       loczsum += local_coords.Z();
1218       // add the association between this cluster key and this hitkey to the table
1219     }  // mapiter
1220     if (mClusHitsVerbose)
1221     {
1222       if (Verbosity() > 10)
1223       {
1224         for (auto& hit : m_iphi)
1225         {
1226           std::cout << " m_phi(" << hit.first << " : " << hit.second << ") " << std::endl;
1227         }
1228       }
1229       for (auto& hit : m_iphi)
1230       {
1231         mClusHitsVerbose->addPhiHit(hit.first, (float) hit.second);
1232       }
1233       for (auto& hit : m_it)
1234       {
1235         mClusHitsVerbose->addZHit(hit.first, (float) hit.second);
1236       }
1237       for (auto& hit : m_iphiCut)
1238       {
1239         mClusHitsVerbose->addPhiCutHit(hit.first, (float) hit.second);
1240       }
1241       for (auto& hit : m_itCut)
1242       {
1243         mClusHitsVerbose->addZCutHit(hit.first, (float) hit.second);
1244       }
1245     }
1246 
1247     // This is the local position
1248     locclusx = locxsum / npixels;
1249     locclusz = loczsum / npixels;
1250 
1251     const double pitch = layergeom->get_pixel_x();
1252     const double length = layergeom->get_pixel_z();
1253     const double phisize = phibins.size() * pitch;
1254     const double zsize = zbins.size() * length;
1255 
1256     /* if (true) { */
1257     if (Verbosity() > 0)
1258     {
1259       std::cout << " MvtxClusterizer: cluskey " << ckey << " layer " << layer << " rad " << layergeom->get_radius() << " phibins " << phibins.size() << " pitch " << pitch << " phisize " << phisize
1260                 << " zbins " << zbins.size() << " length " << length << " zsize " << zsize
1261                 << " local x " << locclusx << " local y " << locclusz
1262                 << std::endl;
1263     }
1264 
1265     // ok force it use use cluster version v4 for now (Valgrind is not happy with application of v5)
1266     /* if (m_cluster_version==4){ */
1267     if (m_cluster_version == 4)
1268     {
1269       auto clus = std::make_unique<TrkrClusterv4>();
1270       clus->setAdc(npixels);
1271       clus->setLocalX(locclusx);
1272       clus->setLocalY(locclusz);
1273 
1274       clus->setPhiSize(phibins.size());
1275       clus->setZSize(zbins.size());
1276       // All silicon surfaces have a 1-1 map to hitsetkey.
1277       // So set subsurface key to 0
1278       clus->setSubSurfKey(0);
1279 
1280       if (Verbosity() > 2)
1281       {
1282         clus->identify();
1283       }
1284 
1285       // get the count of how many clusters have allready been added to this hitsetkey (possibly from other embedded tracks tracks)
1286       m_hitsetkey_cnt.try_emplace(hitsetkey, 0);
1287       unsigned int& cnt = m_hitsetkey_cnt[hitsetkey];
1288       ckey = TrkrDefs::genClusKey(hitsetkey, cnt);
1289       m_truthclusters->addClusterSpecifyKey(ckey, clus.release());
1290       m_current_track->addCluster(ckey);
1291       if (mClusHitsVerbose)
1292       {
1293         mClusHitsVerbose->push_hits(ckey);
1294       }
1295       ++cnt;
1296     }
1297     else
1298     {
1299       std::cout << PHWHERE << std::endl;
1300       std::cout << "Error: only cluster version 4 allowed." << std::endl;
1301     }  // loop over hitsets
1302   }
1303   m_truth_hits->Reset();
1304   prior_g4hit = nullptr;
1305   return;
1306 }
1307 
1308 void PHG4MvtxHitReco::makePixelMask(hitMask& aMask, const std::string& dbName, const std::string& totalPixelsToMask)
1309 {
1310   std::string database = CDBInterface::instance()->getUrl(dbName);
1311   CDBTTree* cdbttree = new CDBTTree(database);
1312 
1313   int NPixel = -1;
1314   NPixel = cdbttree->GetSingleIntValue(totalPixelsToMask);
1315 
1316   for (int i = 0; i < NPixel; i++)
1317   {
1318     int Layer = cdbttree->GetIntValue(i, "layer");
1319     int Stave = cdbttree->GetIntValue(i, "stave");
1320     int Chip = cdbttree->GetIntValue(i, "chip");
1321     int Col = cdbttree->GetIntValue(i, "col");
1322     int Row = cdbttree->GetIntValue(i, "row");
1323     if (Verbosity() > VERBOSITY_A_LOT)
1324     {
1325       std::cout << dbName << ": Will mask layer: " << Layer << ", stave: " << Stave << ", chip: " << Chip << ", row: " << Row << ", col: " << Col << std::endl;
1326     }
1327 
1328     TrkrDefs::hitsetkey DeadPixelHitKey = MvtxDefs::genHitSetKey(Layer, Stave, Chip, 0);
1329     TrkrDefs::hitkey DeadHitKey = MvtxDefs::genHitKey(Col, Row);
1330     aMask.push_back({std::make_pair(DeadPixelHitKey, DeadHitKey)});
1331   }
1332 
1333   delete cdbttree;
1334 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!