Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-06 08:18:41

0001 #include "TrackerEventDisplay.h"
0002 
0003 #include <trackbase/ActsGeometry.h>
0004 #include <trackbase/TrackFitUtils.h>
0005 #include <trackbase/TrkrCluster.h>
0006 #include <trackbase/TrkrDefs.h>
0007 #include <trackbase/TrkrHit.h>
0008 #include <trackbase/TrkrHitSetContainer.h>
0009 
0010 #include <trackbase_historic/TrackSeedContainer_v1.h>
0011 
0012 #include <trackbase/TpcDefs.h>
0013 
0014 #include <trackbase/TrkrClusterContainer.h>
0015 #include <trackbase/TrkrClusterHitAssoc.h>
0016 #include <trackbase/TrkrClusterIterationMapv1.h>
0017 #include <trackbase/TrkrHitSet.h>
0018 
0019 #include <g4detectors/PHG4TpcCylinderGeom.h>
0020 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0021 
0022 #include <fun4all/Fun4AllReturnCodes.h>
0023 #include <fun4all/SubsysReco.h>
0024 
0025 #include <phool/getClass.h>
0026 #include <phool/phool.h>
0027 #include <phool/recoConsts.h>
0028 
0029 #include <TVector3.h>
0030 #include <boost/format.hpp>
0031 #include <boost/math/special_functions/sign.hpp>
0032 
0033 #include <cmath>
0034 #include <iomanip>
0035 #include <iostream>
0036 #include <iterator>
0037 #include <limits>
0038 #include <map>
0039 #include <memory>  // for shared_ptr
0040 #include <set>     // for _Rb_tree_cons...
0041 #include <utility>
0042 #include <vector>
0043 
0044 TrackerEventDisplay::TrackerEventDisplay(const std::string& /*name*/, const std::string& filename, const std::string& runnumber, const std::string& date)
0045   : SubsysReco("TrackerEventDisplay")
0046   , _filename(filename)
0047   , _runnumber(runnumber)
0048   , _date(date)
0049 {
0050 }
0051 
0052 int TrackerEventDisplay::Init(PHCompositeNode* /*topNode*/)
0053 {
0054   _ievent = 0;
0055 
0056   return Fun4AllReturnCodes::EVENT_OK;
0057 }
0058 int TrackerEventDisplay::InitRun(PHCompositeNode *topNode)
0059 {
0060   auto geom =
0061       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0062   if (!geom)
0063   {
0064     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0065     return Fun4AllReturnCodes::ABORTRUN;
0066   }
0067   AdcClockPeriod = geom->GetFirstLayerCellGeom()->get_zstep();
0068 
0069   return Fun4AllReturnCodes::EVENT_OK;
0070 }
0071 int TrackerEventDisplay::process_event(PHCompositeNode* topNode)
0072 {
0073   makeJsonFile(topNode);
0074   ++_ievent;
0075   return Fun4AllReturnCodes::EVENT_OK;
0076 }
0077 
0078 int TrackerEventDisplay::End(PHCompositeNode* /*topNode*/)
0079 {
0080   if (Verbosity() > 1)
0081   {
0082     std::cout << "========================= TrackerEventDisplay::End() ============================" << std::endl;
0083     std::cout << " " << _ievent << " events of output written to: " << _filename << std::endl;
0084     std::cout << "===========================================================================" << std::endl;
0085   }
0086 
0087   return Fun4AllReturnCodes::EVENT_OK;
0088 }
0089 
0090 void TrackerEventDisplay::makeJsonFile(PHCompositeNode* topNode)
0091 {
0092   if (Verbosity() > 1)
0093   {
0094     std::cout << "TrackerEventDisplay::makeJsonFile() entered" << std::endl;
0095   }
0096 
0097   ActsGeometry* tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0098   if (!tgeometry)
0099   {
0100     std::cout << "No Acts geometry on node tree. Can't  continue."
0101               << std::endl;
0102     return;
0103   }
0104 
0105   PHG4TpcCylinderGeomContainer* geom_container =
0106       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0107   if (!geom_container)
0108   {
0109     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0110     return;
0111   }
0112 
0113   //--------------------
0114   // fill the Hit json
0115   //--------------------
0116 
0117   if (_hit)
0118   {
0119     bool firstHit = true;
0120 
0121     outdata.open((_filename + "_event" + std::to_string(_ievent) + "_hits.json").c_str(), std::ofstream::out | std::ofstream::trunc);
0122 
0123     if (!outdata)
0124     {  // file couldn't be opened
0125       std::cout << "ERROR: file could not be opened" << std::endl;
0126       exit(1);
0127     }
0128 
0129     outdata << "{\n    \"EVENT\": {\n        \"runid\":" << _runnumber << ", \n        \"evtid\": 1, \n        \"time\": 0, \n \"timeStr\": \"2023-08-23, 15:23:30 EST\", \n       \"type\": \"Cosmics\", \n        \"s_nn\": 0, \n        \"B\": 0.0,\n        \"pv\": [0,0,0],\n  \"runstats\": [ \n  \"sPHENIX Time Projection Chamber\", \"" << _date << ", Run " << _runnumber << " - Event " << _ievent << "\", \"All Hits in Event\"] \n   },\n"
0130             << std::endl;
0131 
0132     outdata << "    \"META\": {\n       \"HITS\": {\n          \"INNERTRACKER\": {\n              \"type\": \"3D\",\n              \"options\": {\n              \"size\": 2,\n              \"color\": 16777215\n              } \n          },\n"
0133             << std::endl;
0134     outdata << "          \"TRACKHITS\": {\n              \"type\": \"3D\",\n              \"options\": {\n              \"size\": 2,\n              \"transparent\": 0.5,\n              \"color\": 16777215\n              } \n          },\n"
0135             << std::endl;
0136     outdata << "    \"JETS\": {\n        \"type\": \"JET\",\n        \"options\": {\n            \"rmin\": 0,\n            \"rmax\": 78,\n            \"emin\": 0,\n            \"emax\": 30,\n            \"color\": 16777215,\n            \"transparent\": 0.5 \n        }\n    }\n        }\n    }\n," << std::endl;
0137     outdata << "    \"HITS\": {\n        \"CEMC\":[{\"eta\": 0, \"phi\": 0, \"e\": 0}\n            ],\n        \"HCALIN\": [{\"eta\": 0, \"phi\": 0, \"e\": 0}\n            ],\n        \"HCALOUT\": [{\"eta\": 0, \"phi\": 0, \"e\": 0}\n \n            ],\n\n"
0138             << std::endl;
0139     outdata << "    \"TRACKHITS\": [\n\n ";
0140 
0141     auto m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0142 
0143     if (Verbosity() >= 1)
0144     {
0145       std::cout << "Filling hit json" << std::endl;
0146     }
0147     // need things off of the DST...
0148     TrkrHitSetContainer* hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0149 
0150     if (hitmap)
0151     {
0152       TrkrHitSetContainer::ConstRange all_hitsets = hitmap->getHitSets();
0153       for (TrkrHitSetContainer::ConstIterator iter = all_hitsets.first;
0154            iter != all_hitsets.second;
0155            ++iter)
0156       {
0157         const TrkrDefs::hitsetkey hitset_key = iter->first;
0158         TrkrHitSet* hitset = iter->second;
0159         // get all hits for this hitset
0160         TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0161         for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0162              hitr != hitrangei.second;
0163              ++hitr)
0164         {
0165           TrkrDefs::hitkey hit_key = hitr->first;
0166           TrkrHit* hit = hitr->second;
0167           // float event = _ievent;
0168           // float hitID = hit_key;
0169           // float e = hit->getEnergy();
0170           float adc = hit->getAdc();
0171           float layer_local = TrkrDefs::getLayer(hitset_key);
0172           // float sector = TpcDefs::getSectorId(hitset_key);
0173           float side = TpcDefs::getSide(hitset_key);
0174           // float cellID = 0;
0175           // float ecell = hit->getAdc();
0176 
0177           float phibin = std::numeric_limits<float>::quiet_NaN();
0178           float tbin = std::numeric_limits<float>::quiet_NaN();
0179           // float phi = std::numeric_limits<float>::quiet_NaN();
0180           float phi_center = std::numeric_limits<float>::quiet_NaN();
0181           float x = std::numeric_limits<float>::quiet_NaN();
0182           float y = std::numeric_limits<float>::quiet_NaN();
0183           float z = std::numeric_limits<float>::quiet_NaN();
0184 
0185           if (TrkrDefs::getTrkrId(hitset_key) == TrkrDefs::TrkrId::tpcId)
0186           {
0187             PHG4TpcCylinderGeom* GeoLayer_local = geom_container->GetLayerCellGeom(layer_local);
0188             double radius = GeoLayer_local->get_radius();
0189             phibin = (float) TpcDefs::getPad(hit_key);
0190             tbin = (float) TpcDefs::getTBin(hit_key);
0191             // phi = GeoLayer_local->get_phicenter(phibin);
0192 
0193             double zdriftlength = tbin * m_tGeometry->get_drift_velocity() * AdcClockPeriod;
0194             // convert z drift length to z position in the TPC
0195             //      std::cout << " tbin: " << tbin << " vdrift " <<m_tGeometry->get_drift_velocity() << " l drift: " << zdriftlength  <<std::endl;
0196             unsigned short NTBins = (unsigned short) GeoLayer_local->get_zbins();
0197             double m_tdriftmax = AdcClockPeriod * NTBins / 2.0;
0198             double clusz = (m_tdriftmax * m_tGeometry->get_drift_velocity()) - zdriftlength;
0199             if (side == 0)
0200             {
0201               clusz = -clusz;
0202             }
0203             z = clusz;
0204             phi_center = GeoLayer_local->get_phicenter(phibin, side);
0205             x = radius * cos(phi_center);
0206             y = radius * sin(phi_center);
0207 
0208             std::stringstream spts;
0209 
0210             if (firstHit)
0211             {
0212               firstHit = false;
0213             }
0214             else
0215             {
0216               spts << ",";
0217             }
0218 
0219             spts << "{ \"x\": ";
0220             spts << x;
0221             spts << ", \"y\": ";
0222             spts << y;
0223             spts << ", \"z\": ";
0224             spts << z;
0225             spts << ", \"e\": ";
0226             spts << adc;
0227             spts << "}";
0228 
0229             outdata << (boost::format("%1%") % spts.str());
0230             spts.clear();
0231             spts.str("");
0232           }
0233         }
0234       }
0235     }
0236     outdata << "],\n    \"JETS\": [\n         ]\n    }," << std::endl;
0237     outdata << "\"TRACKS\": {" << std::endl;
0238     outdata << "\""
0239             << "INNERTRACKER"
0240             << "\": [";
0241     outdata << "]" << std::endl;
0242     outdata << "}" << std::endl;
0243     outdata << "}" << std::endl;
0244     outdata.close();
0245   }
0246 
0247   //------------------------
0248   // fill the Cluster JSON
0249   //------------------------
0250 
0251   if (_cluster)
0252   {
0253     bool firstHit = true;
0254 
0255     outdata.open((_filename + "_event" + std::to_string(_ievent) + "_clusters.json").c_str(), std::ofstream::out | std::ofstream::trunc);
0256 
0257     if (!outdata)
0258     {  // file couldn't be opened
0259       std::cout << "ERROR: file could not be opened" << std::endl;
0260       exit(1);
0261     }
0262 
0263     outdata << "{\n    \"EVENT\": {\n        \"runid\":" << _runnumber << ", \n        \"evtid\": 1, \n        \"time\": 0, \n \"timeStr\": \"2023-08-23, 15:23:30 EST\", \n       \"type\": \"Cosmics\", \n        \"s_nn\": 0, \n        \"B\": 0.0,\n        \"pv\": [0,0,0],\n  \"runstats\": [ \n  \"sPHENIX Time Projection Chamber\", \"" << _date << ", Run " << _runnumber << " - Event " << _ievent << "\", \"All Clusters in Event\"] \n   },\n"
0264             << std::endl;
0265 
0266     outdata << "    \"META\": {\n       \"HITS\": {\n          \"INNERTRACKER\": {\n              \"type\": \"3D\",\n              \"options\": {\n              \"size\": 2,\n              \"color\": 16777215\n              } \n          },\n"
0267             << std::endl;
0268     outdata << "          \"TRACKHITS\": {\n              \"type\": \"3D\",\n              \"options\": {\n              \"size\": 2,\n              \"transparent\": 0.5,\n              \"color\": 16777215\n              } \n          },\n"
0269             << std::endl;
0270     outdata << "    \"JETS\": {\n        \"type\": \"JET\",\n        \"options\": {\n            \"rmin\": 0,\n            \"rmax\": 78,\n            \"emin\": 0,\n            \"emax\": 30,\n            \"color\": 16777215,\n            \"transparent\": 0.5 \n        }\n    }\n        }\n    }\n," << std::endl;
0271     outdata << "    \"HITS\": {\n        \"CEMC\":[{\"eta\": 0, \"phi\": 0, \"e\": 0}\n            ],\n        \"HCALIN\": [{\"eta\": 0, \"phi\": 0, \"e\": 0}\n            ],\n        \"HCALOUT\": [{\"eta\": 0, \"phi\": 0, \"e\": 0}\n \n            ],\n\n"
0272             << std::endl;
0273     outdata << "    \"TRACKHITS\": [\n\n ";
0274 
0275     if (Verbosity() > 1)
0276     {
0277       std::cout << "Filling cluster json " << std::endl;
0278     }
0279 
0280     // need things off of the DST...
0281     TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
0282     if (!clustermap)
0283     {
0284       clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0285     }
0286 
0287     if (Verbosity() > 1)
0288     {
0289       if (clustermap != nullptr)
0290       {
0291         std::cout << "got clustermap" << std::endl;
0292       }
0293       else
0294       {
0295         std::cout << "no clustermap" << std::endl;
0296       }
0297     }
0298 
0299     if (clustermap)
0300     {
0301       for (const auto& hitsetkey : clustermap->getHitSetKeys())
0302       {
0303         auto range = clustermap->getClusters(hitsetkey);
0304         for (auto iter = range.first; iter != range.second; ++iter)
0305         {
0306           TrkrDefs::cluskey cluster_key = iter->first;
0307           TrkrCluster* cluster = clustermap->findCluster(cluster_key);
0308 
0309           Acts::Vector3 cglob;
0310           cglob = tgeometry->getGlobalPosition(cluster_key, cluster);
0311           float x = cglob(0);
0312           float y = cglob(1);
0313           float z = cglob(2);
0314           float adc = cluster->getAdc();
0315 
0316           std::stringstream spts;
0317 
0318           if (firstHit)
0319           {
0320             firstHit = false;
0321           }
0322           else
0323           {
0324             spts << ",";
0325           }
0326 
0327           spts << "{ \"x\": ";
0328           spts << x;
0329           spts << ", \"y\": ";
0330           spts << y;
0331           spts << ", \"z\": ";
0332           spts << z;
0333           spts << ", \"e\": ";
0334           spts << adc;
0335           spts << "}";
0336 
0337           outdata << (boost::format("%1%") % spts.str());
0338           spts.clear();
0339           spts.str("");
0340         }
0341       }
0342     }
0343     outdata << "],\n    \"JETS\": [\n         ]\n    }," << std::endl;
0344     outdata << "\"TRACKS\": {" << std::endl;
0345     outdata << "\""
0346             << "INNERTRACKER"
0347             << "\": [";
0348     outdata << "]" << std::endl;
0349     outdata << "}" << std::endl;
0350     outdata << "}" << std::endl;
0351     outdata.close();
0352   }
0353   return;
0354 }