Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TPCHitTrackDisplay.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TPCHitTrackDisplay.cc
1 #include "TPCHitTrackDisplay.h"
2 
3 #include <fun4all/Fun4AllReturnCodes.h>
4 #include <fun4all/PHTFileServer.h>
5 #include <fun4all/Fun4AllServer.h>
6 
7 #include <centrality/CentralityInfo.h>
8 #include <centrality/CentralityInfov1.h>
9 
10 #include <g4main/PHG4HitContainer.h>
11 #include <g4main/PHG4Hit.h>
12 #include <g4main/PHG4TruthInfoContainer.h>
13 #include <g4main/PHG4VtxPoint.h>
14 #include <g4main/PHG4Particle.h>
15 
16 #include <phool/PHCompositeNode.h>
17 #include <phool/PHIODataNode.h>
18 #include <phool/PHNode.h>
19 #include <phool/PHNodeIterator.h>
20 #include <phool/PHObject.h>
21 #include <phool/getClass.h>
22 #include <phool/phool.h>
23 
24 #include <trackbase/TrkrCluster.h>
25 
26 #include <phgeom/PHGeomUtility.h>
27 
28 #include <TTree.h>
29 #include <TH2D.h>
30 #include <TVector3.h>
31 #include <TRandom3.h>
32 #include <TMath.h>
33 
34 #include <trackbase_historic/SvtxTrackMap.h>
35 #include <jetbase/Jet.h>
36 #include <jetbase/JetContainer.h>
37 #include <jetbase/JetContainerv1.h>
38 #include <g4main/PHG4Utils.h>
39 
40 #include <iostream>
41 #include <cassert>
42 #include <fstream>
43 #include <sstream>
44 #include <boost/format.hpp>
45 #include <boost/math/special_functions/sign.hpp>
46 
47 /*************************************************************/
48 /* TPC Event Display Generator */
49 /* Thomas Marshall,Aditya Dash,Ejiro Umaka */
50 /*rosstom@ucla.edu,aditya55@physics.ucla.edu,eumaka1@bnl.gov */
51 /*************************************************************/
52 
53 
54 using namespace std;
55 
56 //----------------------------------------------------------------------------//
57 //-- Constructor:
58 //-- simple initialization
59 //----------------------------------------------------------------------------//
60 
61 TPCHitTrackDisplay::TPCHitTrackDisplay(const string &name) :
62  SubsysReco(name)
63  , m_cut_ADC(75.0)
64  , m_trackless_clusters(true)
65  , _pdgcode(0)
66  , _fileName(name)
67 {
68  //initialize
69  _event = 0;
70 }
71 
72 //----------------------------------------------------------------------------//
73 //-- Init():
74 //-- Intialize all histograms, trees, and ntuples
75 //----------------------------------------------------------------------------//
76 int TPCHitTrackDisplay::Init(PHCompositeNode *topNode) {
77 
78  return Fun4AllReturnCodes::EVENT_OK;
79 }
80 
81 int TPCHitTrackDisplay::InitRun(PHCompositeNode *topNode)
82 {
83  return Fun4AllReturnCodes::EVENT_OK;
84 }
85 
86 //----------------------------------------------------------------------------//
87 //-- process_event():
88 //-- Call user instructions for every event.
89 //-- This function contains the analysis structure.
90 //----------------------------------------------------------------------------//
91 
92 int TPCHitTrackDisplay::process_event(PHCompositeNode *topNode)
93 {
94  _event++;
95  SimulationOut(topNode);
96 
97  return Fun4AllReturnCodes::EVENT_OK;
98 }
99 
100 //----------------------------------------------------------------------------//
101 //-- End():
102 //-- End method, wrap everything up
103 //----------------------------------------------------------------------------//
104 
105 int TPCHitTrackDisplay::EndRun(const int runnumber)
106 {
107 
108  return Fun4AllReturnCodes::EVENT_OK;
109 
110 }
111 
112 //Produce json file from raw TPC csv data
113 // Will adjust this soon, once TrkrHitSetContainer is updated
114 /*
115 void TPCHitTrackDisplay::TPCRawOut(PHCompositeNode *topNode)
116 {
117 
118  cout << _event << endl;
119 
120  outfile1 << "\"TRACKS\": {" << endl;
121  outfile1 <<"\""<<"INNERTRACKER"<<"\": [";
122 
123  bool first = true;
124 
125  stringstream spts;
126  float c = NAN;
127 
128  float x = 0.; float y = 0.; float z = 0.;
129  float px = 0.; float py = 0.; float pz = 0.;
130  TVector3 pos; TVector3 mom;
131 
132  for
133 
134  for (SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end(); ++iter)
135  {
136  SvtxTrack* track = iter->second;
137  px = track->get_px();
138  py = track->get_py();
139  pz = track->get_pz();
140  mom.SetX(px); mom.SetY(py); mom.SetZ(pz);
141  c = track->get_charge();
142 
143 
144  std::vector<TrkrDefs::cluskey> clusters;
145  auto siseed = track->get_silicon_seed();
146  if(siseed)
147  {
148  for (auto iter = siseed->begin_cluster_keys(); iter != siseed->end_cluster_keys(); ++iter)
149  {
150  TrkrDefs::cluskey cluster_key = *iter;
151  clusters.push_back(cluster_key);
152  }
153  }
154 
155  auto tpcseed = track->get_tpc_seed();
156  if(tpcseed)
157  {
158  for (auto iter = tpcseed->begin_cluster_keys(); iter != tpcseed->end_cluster_keys(); ++iter)
159  {
160  TrkrDefs::cluskey cluster_key = *iter;
161  clusters.push_back(cluster_key);
162  }
163  }
164 
165 
166  for(unsigned int iclus = 0; iclus < clusters.size(); ++iclus)
167  {
168  TrkrDefs::cluskey cluster_key = clusters[iclus];
169  TrkrCluster* cluster = clusterContainer->findCluster(cluster_key);
170  if(!cluster) continue;
171 
172  Acts::Vector3 globalClusterPosition = geometry->getGlobalPosition(cluster_key, cluster);
173  x = globalClusterPosition(0);
174  y = globalClusterPosition(1);
175  z = globalClusterPosition(2);
176  pos.SetX(x); pos.SetY(y); pos.SetZ(z);
177 
178  if (first)
179  {
180  first = false;
181  }
182 
183  else
184  spts << ",";
185  spts << "[";
186  spts << pos.x();
187  spts << ",";
188  spts << pos.y();
189  spts << ",";
190  spts << pos.z();
191  spts << "]";
192 
193  first = true;
194  }
195 
196  outfile1
197  << (boost::format(
198  "{ \"pt\": %1%, \"e\": %2%, \"p\": %3%, \"c\": %4%, \"pdgcode \": %5%, \"pts\":[ %6% ]},")
199  % mom.Pt() % mom.PseudoRapidity() % mom.Phi() % c % _pdgcode % spts.str() ) << endl;
200  spts.clear();
201  spts.str("");
202 }
203 
204 
205  outfile1 << "]" << endl;
206  outfile1 << "}" << endl;
207 
208 
209 
210 return;
211 
212 }
213 */
214 
215 // write json file from simulation data for silicon and tpc clusters/tracks
216 void TPCHitTrackDisplay::SimulationOut(PHCompositeNode *topNode)
217 {
218  stringstream fname;
219  fname << _fileName << "event" << _event << ".json";
220  fstream outfile1(fname.str(), ios_base::out);
221 
222  cout << _event << endl;
223 
224  TrkrClusterContainer *clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
225 
226  CentralityInfov1 *cent = findNode::getClass<CentralityInfov1>(topNode, "CentralityInfo");
227  if (!cent)
228  {
229  std::cout << " ERROR -- can't find CentralityInfo node" << std::endl;
230  return;
231  }
232 
233  int cent_index = cent->get_centile(CentralityInfo::PROP::bimp);
234 
235 
236  SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
237 
238  ActsGeometry *geometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
239  if(!geometry)
240  {
241  std::cout << PHWHERE << "No Acts geometry on node tree. Can't continue."
242  << std::endl;
243  }
244 
245  std::cout<< "This event has centrality: \t" << cent_index << std::endl;
246 
247 
248  // Header information for the json file
249 
250  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" << endl;
251 
252  outfile1 << " \"META\": {\n \"HITS\": {\n \"INNERTRACKER\": {\n \"type\": \"3D\",\n \"options\": {\n \"size\": 5,\n \"color\": 16777215\n } \n },\n" << endl;
253  outfile1 << " \"TRACKHITS\": {\n \"type\": \"3D\",\n \"options\": {\n \"size\": 5,\n \"transparent\": 0.5,\n \"color\": 16777215\n } \n },\n" << endl;
254  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," << endl;
255  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" << endl;
256  outfile1 << " \"TRACKHITS\": [\n\n ";
257 
258  std::vector<TrkrDefs::cluskey> usedClusters;
259 
260  // Fill usedClusters with cluster keys that are associated with a reconstructed track
261  for (SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end(); ++iter)
262  {
263  if (!m_trackless_clusters) break;
264  SvtxTrack* track = iter->second;
265  std::vector<TrkrDefs::cluskey> clusters;
266  auto siseed = track->get_silicon_seed();
267  if(siseed)
268  {
269  for (auto iter = siseed->begin_cluster_keys(); iter != siseed->end_cluster_keys(); ++iter)
270  {
271  TrkrDefs::cluskey cluster_key = *iter;
272  usedClusters.push_back(cluster_key);
273  }
274  }
275 
276  auto tpcseed = track->get_tpc_seed();
277  if(tpcseed)
278  {
279  for (auto iter = tpcseed->begin_cluster_keys(); iter != tpcseed->end_cluster_keys(); ++iter)
280  {
281  TrkrDefs::cluskey cluster_key = *iter;
282  usedClusters.push_back(cluster_key);
283  }
284  }
285  }
286 
287  bool firstClus = true;
288  int counter = 0;
289  stringstream spts;
290  float c = NAN;
291 
292  float x = 0.; float y = 0.; float z = 0.;
293  float px = 0.; float py = 0.; float pz = 0.;
294  TVector3 pos; TVector3 mom;
295 
296  // iterate over all clusters and write out only those without an associated track
297  for(const auto& hitsetkey:clusterContainer->getHitSetKeys())
298  {
299  if (!m_trackless_clusters) break; // cut on whether or not to display trackless clusters
300  auto range = clusterContainer->getClusters(hitsetkey);
301  for( auto iter = range.first; iter != range.second; ++iter )
302  {
303  counter = counter + 1;
304  TrkrDefs::cluskey cluster_key = iter->first;
305  bool clusFromTrack = false;
306  for(unsigned int iclus = 0; iclus < usedClusters.size(); ++iclus)
307  {
308  TrkrDefs::cluskey trackClusKey = usedClusters[iclus];
309  if (trackClusKey == cluster_key)
310  {
311  clusFromTrack = true;
312  break;
313  }
314  }
315  if (clusFromTrack) continue;
316  TrkrCluster* cluster = iter->second;
317  if(!cluster) continue;
318 
319  unsigned int clusADC = cluster->getAdc();
320  if (clusADC < m_cut_ADC) continue; // ADC check on cluster to not over-saturate display
321 
322  Acts::Vector3 globalClusterPosition = geometry->getGlobalPosition(cluster_key, cluster);
323  x = globalClusterPosition(0);
324  y = globalClusterPosition(1);
325  z = globalClusterPosition(2);
326  pos.SetX(x); pos.SetY(y); pos.SetZ(z);
327 
328  if (firstClus) firstClus = false;
329  else spts << ",";
330 
331  spts << "{ \"x\": ";
332  spts << pos.x();
333  spts << ", \"y\": ";
334  spts << pos.y();
335  spts << ", \"z\": ";
336  spts << pos.z();
337  spts << ", \"e\": ";
338  spts << clusADC;
339  spts << "}";
340 
341  outfile1 << (boost::format("%1%") % spts.str());
342  spts.clear();
343  spts.str("");
344  }
345  }
346 
347  outfile1 << "],\n \"JETS\": [\n ]\n }," << endl;
348 
349  outfile1 << "\"TRACKS\": {" << endl;
350  outfile1 <<"\""<<"INNERTRACKER"<<"\": [";
351 
352  firstClus = true;
353  bool firstTrack = true;
354 
355  //tracking
356 
357  c = NAN;
358 
359  x = 0.; y = 0.; z = 0.;
360  px = 0.; py = 0.; pz = 0.;
361 
362  // iterate over tracks and write to file all tracks and their associated clusters
363  for (SvtxTrackMap::Iter iter = trackmap->begin(); iter != trackmap->end(); ++iter)
364  {
365  SvtxTrack* track = iter->second;
366  px = track->get_px();
367  py = track->get_py();
368  pz = track->get_pz();
369  mom.SetX(px); mom.SetY(py); mom.SetZ(pz);
370  c = track->get_charge();
371 
372  std::vector<TrkrDefs::cluskey> clusters;
373  auto siseed = track->get_silicon_seed();
374  if(siseed)
375  {
376  for (auto iter = siseed->begin_cluster_keys(); iter != siseed->end_cluster_keys(); ++iter)
377  {
378  TrkrDefs::cluskey cluster_key = *iter;
379  clusters.push_back(cluster_key);
380  }
381  }
382 
383  auto tpcseed = track->get_tpc_seed();
384  if(tpcseed)
385  {
386  for (auto iter = tpcseed->begin_cluster_keys(); iter != tpcseed->end_cluster_keys(); ++iter)
387  {
388  TrkrDefs::cluskey cluster_key = *iter;
389  clusters.push_back(cluster_key);
390  }
391  }
392 
393 
394  for(unsigned int iclus = 0; iclus < clusters.size(); ++iclus)
395  {
396  TrkrDefs::cluskey cluster_key = clusters[iclus];
397  TrkrCluster* cluster = clusterContainer->findCluster(cluster_key);
398  if(!cluster) continue;
399 
400  Acts::Vector3 globalClusterPosition = geometry->getGlobalPosition(cluster_key, cluster);
401  x = globalClusterPosition(0);
402  y = globalClusterPosition(1);
403  z = globalClusterPosition(2);
404  pos.SetX(x); pos.SetY(y); pos.SetZ(z);
405 
406  if (firstClus)
407  {
408  firstClus = false;
409  }
410 
411  else
412  spts << ",";
413  spts << "[";
414  spts << pos.x();
415  spts << ",";
416  spts << pos.y();
417  spts << ",";
418  spts << pos.z();
419  spts << "]";
420 
421  }
422 
423  firstClus = true;
424  if (firstTrack)
425  {
426  outfile1
427  << (boost::format(
428  "{ \"pt\": %1%, \"e\": %2%, \"p\": %3%, \"c\": %4%, \"pdgcode \": %5%, \"centrality \": %6%, \"pts\":[ %7% ]}")
429  % mom.Pt() % mom.PseudoRapidity() % mom.Phi() % c % _pdgcode % cent_index % spts.str() ) << endl;
430  spts.clear();
431  spts.str("");
432  firstTrack = false;
433  }
434  else
435  {
436  outfile1
437  << (boost::format(
438  ",{ \"pt\": %1%, \"e\": %2%, \"p\": %3%, \"c\": %4%, \"pdgcode \": %5%, \"centrality \": %6%, \"pts\":[ %7% ]}")
439  % mom.Pt() % mom.PseudoRapidity() % mom.Phi() % c % _pdgcode % cent_index % spts.str() ) << endl;
440  spts.clear();
441  spts.str("");
442  }
443  }
444 
445 
446  outfile1 << "]" << endl;
447  outfile1 << "}" << endl;
448  outfile1 << "}" << endl;
449 
450  usedClusters.clear();
451 
452  return;
453 
454 }
455