Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:22:08 

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