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