packages/TPCHitTrackDisplay/TPCHitTrackDisplay.cc

0001 #include "TPCHitTrackDisplay.h"
0002
0003 #include <centrality/CentralityInfo.h>
0004
0005 #include <trackbase/ActsGeometry.h>
0006 #include <trackbase/TrkrCluster.h>
0007 #include <trackbase/TrkrClusterContainer.h>
0008 #include <trackbase/TrkrDefs.h>
0009
0010 #include <trackbase_historic/SvtxTrack.h>
0011 #include <trackbase_historic/SvtxTrackMap.h>
0012 #include <trackbase_historic/TrackSeed.h>
0013
0014 #include <phgeom/PHGeomUtility.h>
0015
0016 #include <fun4all/Fun4AllReturnCodes.h>
0017
0018 #include <phool/getClass.h>
0019 #include <phool/phool.h>
0020
0021 #include <TVector3.h>
0022
0023 #include <format>
0024 #include <fstream>
0025 #include <iostream>
0026 #include <limits>
0027 #include <sstream>
0028
0029 /*************************************************************/
0030 /*              TPC Event Display Generator                  */
0031 /*         Thomas Marshall,Aditya Dash,Ejiro Umaka           */
0032 /*rosstom@ucla.edu,aditya55@physics.ucla.edu,eumaka1@bnl.gov */
0033 /*************************************************************/
0034
0035 //----------------------------------------------------------------------------//
0036 //-- Constructor:
0037 //--  simple initialization
0038 //----------------------------------------------------------------------------//
0039
0040 TPCHitTrackDisplay::TPCHitTrackDisplay(const std::string& name)
0041   : SubsysReco(name)
0042   , m_cut_ADC(75.0)
0043   , m_trackless_clusters(true)
0044   , _fileName(name)
0045 {
0046 }
0047
0048 //----------------------------------------------------------------------------//
0049 //-- process_event():
0050 //--   Call user instructions for every event.
0051 //--   This function contains the analysis structure.
0052 //----------------------------------------------------------------------------//
0053
0054 int TPCHitTrackDisplay::process_event(PHCompositeNode* topNode)
0055 {
0056   _event++;
0057   SimulationOut(topNode);
0058
0059   return Fun4AllReturnCodes::EVENT_OK;
0060 }
0061
0062 // Produce json file from raw TPC csv data
0063 //  Will adjust this soon, once TrkrHitSetContainer is updated
0064 /*
0065 void TPCHitTrackDisplay::TPCRawOut(PHCompositeNode *topNode)
0066 {
0067
0068     std::cout << _event << std::endl;
0069
0070     outfile1 << "\"TRACKS\": {" << std::endl;
0071     outfile1 <<"\""<<"INNERTRACKER"<<"\": [";
0072
0073     bool first = true;
0074
0075     std::stringstream spts;
0076     float c = std::numeric_limits<float>::quiet_NaN();
0077
0078     float x = 0.; float y = 0.; float z = 0.;
0079     float px = 0.; float py = 0.; float pz = 0.;
0080     TVector3 pos; TVector3 mom;
0081
0082     for
0083
0084     for (SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end(); ++iter)
0085     {
0086       SvtxTrack* track = iter->second;
0087       px = track->get_px();
0088       py = track->get_py();
0089       pz = track->get_pz();
0090       mom.SetX(px); mom.SetY(py); mom.SetZ(pz);
0091       c = track->get_charge();
0092
0093
0094       std::vector<TrkrDefs::cluskey> clusters;
0095       auto siseed = track->get_silicon_seed();
0096       if(siseed)
0097       {
0098         for (auto iter = siseed->begin_cluster_keys(); iter != siseed->end_cluster_keys(); ++iter)
0099             {
0100                   TrkrDefs::cluskey cluster_key = *iter;
0101                   clusters.push_back(cluster_key);
0102             }
0103         }
0104
0105       auto tpcseed = track->get_tpc_seed();
0106       if(tpcseed)
0107       {
0108         for (auto iter = tpcseed->begin_cluster_keys(); iter != tpcseed->end_cluster_keys(); ++iter)
0109             {
0110                   TrkrDefs::cluskey cluster_key = *iter;
0111                   clusters.push_back(cluster_key);
0112             }
0113       }
0114
0115
0116      for(unsigned int iclus = 0; iclus < clusters.size(); ++iclus)
0117       {
0118         TrkrDefs::cluskey cluster_key = clusters[iclus];
0119         TrkrCluster* cluster = clusterContainer->findCluster(cluster_key);
0120         if(!cluster) continue;
0121
0122         Acts::Vector3 globalClusterPosition = geometry->getGlobalPosition(cluster_key, cluster);
0123         x = globalClusterPosition(0);
0124         y = globalClusterPosition(1);
0125         z = globalClusterPosition(2);
0126         pos.SetX(x); pos.SetY(y); pos.SetZ(z);
0127
0128         if (first)
0129         {
0130             first = false;
0131         }
0132
0133         else
0134         spts << ",";
0135         spts << "[";
0136         spts << pos.x();
0137         spts << ",";
0138         spts << pos.y();
0139         spts << ",";
0140         spts << pos.z();
0141         spts << "]";
0142
0143         first = true;
0144     }
0145
0146      outfile1
0147         << std::format("{ \"pt\": {}, \"e\": {}, \"p\": {}, \"c\": {}, \"pdgcode \": {}, \"pts\":[ {} ]},",
0148           mom.Pt(), mom.PseudoRapidity(), mom.Phi(), c, _pdgcode, spts.str() ) << std::endl;
0149            spts.str("");
0150 }
0151
0152
0153     outfile1 << "]" << std::endl;
0154     outfile1 << "}" << std::endl;
0155
0156
0157
0158 return;
0159
0160 }
0161 */
0162
0163 // write json file from simulation data for silicon and tpc clusters/tracks
0164 void TPCHitTrackDisplay::SimulationOut(PHCompositeNode* topNode)
0165 {
0166   std::string fname = _fileName + "event" + std::to_string(_event) + ".json";
0167   std::fstream outfile1(fname, std::ios_base::out);
0168
0169   std::cout << _event << std::endl;
0170
0171   TrkrClusterContainer* clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0172
0173   CentralityInfo* cent = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0174   if (!cent)
0175   {
0176     std::cout << " ERROR -- can't find CentralityInfo node" << std::endl;
0177     return;
0178   }
0179
0180   int cent_index = cent->get_centile(CentralityInfo::PROP::bimp);
0181
0182   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0183
0184   ActsGeometry* geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0185   if (!geometry)
0186   {
0187     std::cout << PHWHERE << "No Acts geometry on node tree. Can't  continue."
0188               << std::endl;
0189   }
0190
0191   std::cout << "This event has centrality: \t" << cent_index << std::endl;
0192
0193   // Header information for the json file
0194
0195   outfile1 << "{\n    \"EVENT\": {\n        \"runid\": 1, \n        \"evtid\": 1, \n        \"time\": 0, \n        \"type\": \"Collision\", \n        \"s_nn\": 0, \n        \"B\": 3.0,\n        \"pv\": [0,0,0]  \n    },\n"
0196            << std::endl;
0197
0198   outfile1 << "    \"META\": {\n       \"HITS\": {\n          \"INNERTRACKER\": {\n              \"type\": \"3D\",\n              \"options\": {\n              \"size\": 5,\n              \"color\": 16777215\n              } \n          },\n"
0199            << std::endl;
0200   outfile1 << "          \"TRACKHITS\": {\n              \"type\": \"3D\",\n              \"options\": {\n              \"size\": 5,\n              \"transparent\": 0.5,\n              \"color\": 16777215\n              } \n          },\n"
0201            << std::endl;
0202   outfile1 << "    \"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;
0203   outfile1 << "    \"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"
0204            << std::endl;
0205   outfile1 << "    \"TRACKHITS\": [\n\n ";
0206
0207   std::vector<TrkrDefs::cluskey> usedClusters;
0208
0209   // Fill usedClusters with cluster keys that are associated with a reconstructed track
0210   for (auto& iter : *trackmap)
0211   {
0212     if (!m_trackless_clusters)
0213     {
0214       break;
0215     }
0216     SvtxTrack* track = iter.second;
0217     std::vector<TrkrDefs::cluskey> clusters;
0218     auto *siseed = track->get_silicon_seed();
0219     if (siseed)
0220     {
0221       for (auto iter1 = siseed->begin_cluster_keys(); iter1 != siseed->end_cluster_keys(); ++iter1)
0222       {
0223         TrkrDefs::cluskey cluster_key = *iter1;
0224         usedClusters.push_back(cluster_key);
0225       }
0226     }
0227
0228     auto *tpcseed = track->get_tpc_seed();
0229     if (tpcseed)
0230     {
0231       for (auto iter1 = tpcseed->begin_cluster_keys(); iter1 != tpcseed->end_cluster_keys(); ++iter1)
0232       {
0233         TrkrDefs::cluskey cluster_key = *iter1;
0234         usedClusters.push_back(cluster_key);
0235       }
0236     }
0237   }
0238
0239   bool firstClus = true;
0240   int counter = 0;
0241   std::stringstream spts;
0242   float c = std::numeric_limits<float>::quiet_NaN();
0243
0244   float x = 0.;
0245   float y = 0.;
0246   float z = 0.;
0247   float px = 0.;
0248   float py = 0.;
0249   float pz = 0.;
0250   TVector3 pos;
0251   TVector3 mom;
0252
0253   // iterate over all clusters and write out only those without an associated track
0254   for (const auto& hitsetkey : clusterContainer->getHitSetKeys())
0255   {
0256     if (!m_trackless_clusters)
0257     {
0258       break;  // cut on whether or not to display trackless clusters
0259     }
0260     auto range = clusterContainer->getClusters(hitsetkey);
0261     for (auto iter = range.first; iter != range.second; ++iter)
0262     {
0263       counter = counter + 1;
0264       TrkrDefs::cluskey cluster_key = iter->first;
0265       bool clusFromTrack = false;
0266       for (unsigned long trackClusKey : usedClusters)
0267       {
0268         if (trackClusKey == cluster_key)
0269         {
0270           clusFromTrack = true;
0271           break;
0272         }
0273       }
0274       if (clusFromTrack)
0275       {
0276         continue;
0277       }
0278       TrkrCluster* cluster = iter->second;
0279       if (!cluster)
0280       {
0281         continue;
0282       }
0283
0284       unsigned int clusADC = cluster->getAdc();
0285       if (clusADC < m_cut_ADC)
0286       {
0287         continue;  // ADC check on cluster to not over-saturate display
0288       }
0289
0290       Acts::Vector3 globalClusterPosition = geometry->getGlobalPosition(cluster_key, cluster);
0291       x = globalClusterPosition(0);
0292       y = globalClusterPosition(1);
0293       z = globalClusterPosition(2);
0294       pos.SetX(x);
0295       pos.SetY(y);
0296       pos.SetZ(z);
0297
0298       if (firstClus)
0299       {
0300         firstClus = false;
0301       }
0302       else
0303       {
0304         spts << ",";
0305       }
0306
0307       spts << "{ \"x\": ";
0308       spts << pos.x();
0309       spts << ", \"y\": ";
0310       spts << pos.y();
0311       spts << ", \"z\": ";
0312       spts << pos.z();
0313       spts << ", \"e\": ";
0314       spts << clusADC;
0315       spts << "}";
0316
0317       outfile1 << std::format("{}", spts.str());
0318       spts.str("");
0319     }
0320   }
0321
0322   outfile1 << "],\n    \"JETS\": [\n         ]\n    }," << std::endl;
0323
0324   outfile1 << "\"TRACKS\": {" << std::endl;
0325   outfile1 << "\""
0326            << "INNERTRACKER"
0327            << "\": [";
0328
0329   firstClus = true;
0330   bool firstTrack = true;
0331
0332   // tracking
0333
0334   c = std::numeric_limits<float>::quiet_NaN();
0335
0336   x = 0.;
0337   y = 0.;
0338   z = 0.;
0339   px = 0.;
0340   py = 0.;
0341   pz = 0.;
0342
0343   // iterate over tracks and write to file all tracks and their associated clusters
0344   for (auto& iter : *trackmap)
0345   {
0346     SvtxTrack* track = iter.second;
0347     px = track->get_px();
0348     py = track->get_py();
0349     pz = track->get_pz();
0350     mom.SetX(px);
0351     mom.SetY(py);
0352     mom.SetZ(pz);
0353     c = track->get_charge();
0354
0355     std::vector<TrkrDefs::cluskey> clusters;
0356     auto *siseed = track->get_silicon_seed();
0357     if (siseed)
0358     {
0359       for (auto iter1 = siseed->begin_cluster_keys(); iter1 != siseed->end_cluster_keys(); ++iter1)
0360       {
0361         TrkrDefs::cluskey cluster_key = *iter1;
0362         clusters.push_back(cluster_key);
0363       }
0364     }
0365
0366     auto *tpcseed = track->get_tpc_seed();
0367     if (tpcseed)
0368     {
0369       for (auto iter1 = tpcseed->begin_cluster_keys(); iter1 != tpcseed->end_cluster_keys(); ++iter1)
0370       {
0371         TrkrDefs::cluskey cluster_key = *iter1;
0372         clusters.push_back(cluster_key);
0373       }
0374     }
0375
0376     for (unsigned long cluster_key : clusters)
0377     {
0378       TrkrCluster* cluster = clusterContainer->findCluster(cluster_key);
0379       if (!cluster)
0380       {
0381         continue;
0382       }
0383
0384       Acts::Vector3 globalClusterPosition = geometry->getGlobalPosition(cluster_key, cluster);
0385       x = globalClusterPosition(0);
0386       y = globalClusterPosition(1);
0387       z = globalClusterPosition(2);
0388       pos.SetX(x);
0389       pos.SetY(y);
0390       pos.SetZ(z);
0391
0392       if (firstClus)
0393       {
0394         firstClus = false;
0395       }
0396
0397       else
0398       {
0399         spts << ",";
0400       }
0401       spts << "[";
0402       spts << pos.x();
0403       spts << ",";
0404       spts << pos.y();
0405       spts << ",";
0406       spts << pos.z();
0407       spts << "]";
0408     }
0409
0410     firstClus = true;
0411     if (firstTrack)
0412     {
0413       outfile1 << std::format("{{ \"pt\": {}, \"e\": {}, \"p\": {}, \"c\": {}, \"pdgcode \": {}, \"centrality \": {}, \"pts\": [ {} ] }}",
0414           mom.Pt(), mom.PseudoRapidity(), mom.Phi(), c, _pdgcode, cent_index, spts.str())
0415           << std::endl;
0416       spts.str("");
0417       firstTrack = false;
0418     }
0419     else
0420     {
0421       outfile1 << std::format(",{{ \"pt\": {}, \"e\": {}, \"p\": {}, \"c\": {}, \"pdgcode\": {}, \"centrality\": {}, \"pts\": [ {} ] }}",
0422                   mom.Pt(), mom.PseudoRapidity(), mom.Phi(), c, _pdgcode, cent_index, spts.str())
0423           << std::endl;
0424       spts.str("");
0425     }
0426   }
0427
0428   outfile1 << "]" << std::endl;
0429   outfile1 << "}" << std::endl;
0430   outfile1 << "}" << std::endl;
0431
0432   usedClusters.clear();
0433
0434   return;
0435 }