Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-12-17 09:21:16

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