EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
readerTest.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file readerTest.cc
1 /*******************************************************************************
2  * Copyright (c) The JETSCAPE Collaboration, 2018
3  *
4  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
5  *
6  * For the list of contributors see AUTHORS.
7  *
8  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
9  *
10  * or via email to bugs.jetscape@gmail.com
11  *
12  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
13  * See COPYING for details.
14  ******************************************************************************/
15 // Reader test (focus on graph)
16 
17 #include <iostream>
18 #include <fstream>
19 #include <memory>
20 #include <chrono>
21 #include <thread>
22 
23 #include "gzstream.h"
24 #include "PartonShower.h"
25 #include "JetScapeLogger.h"
26 #include "JetScapeReader.h"
27 #include "JetScapeBanner.h"
28 #include "fjcore.hh"
29 
30 #include <GTL/dfs.h>
31 
32 using namespace std;
33 //using namespace fjcore;
34 
35 using namespace Jetscape;
36 
37 // -------------------------------------
38 
39 // Forward declaration
40 void Show();
42 ostream & operator<<(ostream & ostr, const fjcore::PseudoJet & jet);
43 
44 // -------------------------------------
45 
46 // Create a pdf of the shower graph:
47 // Use with graphviz (on Mac: brew install graphviz --with-app)
48 // in shell: dot GVfile.gv -Tpdf -o outputPDF.pdf
49 // [or you can also use the GraphViz app for Mac Os X (in the "cellar" of homebrew)]
50 
51 // -------------------------------------
52 
53 int main(int argc, char** argv)
54 {
55  JetScapeLogger::Instance()->SetDebug(false);
56  JetScapeLogger::Instance()->SetRemark(false);
57  //SetVerboseLevel (9 a lot of additional debug output ...)
58  //If you want to suppress it: use SetVerboseLevle(0) or max SetVerboseLevle(9) or 10
59  JetScapeLogger::Instance()->SetVerboseLevel(0);
60 
61  cout<<endl;
62  Show();
63 
64  //Do some dummy jetfinding ...
65  fjcore::JetDefinition jet_def(fjcore::antikt_algorithm, 0.7);
66 
67  vector<shared_ptr<PartonShower>> mShowers;
68 
69  //Directly with template: provide the relevant stream
70  //auto reader=make_shared<JetScapeReader<ifstream>>("test_out.dat");
71  //auto reader=make_shared<JetScapeReader<igzstream>>("test_out.dat.gz");
72 
73  // Hide Template (see class declarations in reader/JetScapeReader.h) ...
74  auto reader=make_shared<JetScapeReaderAscii>("test_out.dat");
75  //auto reader=make_shared<JetScapeReaderAsciiGZ>("test_out.dat.gz");
76 
77  // reads in multiple events and multiple shower per event
78  // commentend out so that you get the dot graph file for the first shower in the first event
79  // (add in and the file gets overriden)
80  //while (!reader->Finished())
81  {
82  reader->Next();
83 
84  cout<<"Analyze current event = "<<reader->GetCurrentEvent()<<endl;
85  mShowers=reader->GetPartonShowers();
86 
87  int finals = 0;
88  for (int i=0;i<mShowers.size();i++)
89  {
90  cout<<" Analyze parton shower = "<<i<<endl;
91 
92  mShowers[i]->PrintVertices();
93  mShowers[i]->PrintPartons();
94 
95  finals += mShowers[i]->GetFinalPartonsForFastJet().size();
96 
97  fjcore::ClusterSequence cs(mShowers[i]->GetFinalPartonsForFastJet(), jet_def);
98 
99  vector<fjcore::PseudoJet> jets = fjcore::sorted_by_pt(cs.inclusive_jets(2));
100  cout<<endl;
101  cout<<jet_def.description()<<endl;
102  // Output of found jets ...
103  //cout<<endl;
104  for (int k=0;k<jets.size();k++)
105  cout<<"Anti-kT jet "<<k<<" : "<<jets[k]<<endl;
106  cout<<endl;
107  cout<<"Shower initiating parton : "<<*(mShowers[i]->GetPartonAt(0))<<endl;
108  cout<<endl;
109 
110  AnalyzeGraph(mShowers[i]);
111 
112  if (i==0)
113  {
114  mShowers[i]->SaveAsGV("my_test.gv");
115  mShowers[i]->SaveAsGML("my_test.gml");
116  mShowers[i]->SaveAsGraphML("my_test.graphml");
117  }
118 
119  // wait for 5s
120  //std::this_thread::sleep_for(std::chrono::milliseconds(5000));
121  }
122  cout << " Found " << finals << " final state partons." << endl;
123  auto hadrons = reader->GetHadrons();
124  cout<<"Number of hadrons is: " << hadrons.size() << endl;
125 
126  fjcore::ClusterSequence hcs(reader->GetHadronsForFastJet(), jet_def);
127  vector<fjcore::PseudoJet> hjets = fjcore::sorted_by_pt(hcs.inclusive_jets(2));
128  cout<<"AT HADRONIC LEVEL " << endl;
129  for (int k=0;k<hjets.size();k++)
130  cout<<"Anti-kT jet "<<k<<" : "<<hjets[k]<<endl;
131 
132  // for(unsigned int i=0; i<hadrons.size(); i++) {
133  // cout<<"For Hadron Number "<<i<<" "<< hadrons[i].get()->e() << " "<< hadrons[i].get()->px()<< " "<< hadrons[i].get()->py() << " "<< hadrons[i].get()->pz()<< " "<< hadrons[i].get()->pt()<< endl;
134  // }
135  }
136 
137  reader->Close();
138 }
139 
140 // -------------------------------------
141 
143 {
144  INFO<<"Some GTL graph/shower analysis/dfs search output:";
145 
146  // quick and dirty ...
147  dfs search;
148  search.calc_comp_num(true);
149  search.scan_whole_graph(true);
150  search.start_node();// defaulted to first node ...
151  search.run(*mS);
152 
153  cout<<endl;
154  cout<<"DFS graph search feature from GTL:"<<endl;
155  cout<<"Number of Nodes reached from node 0 = "<< search.number_of_reached_nodes () <<endl;
156  cout<<"Node/Vertex ordering result from DFS:"<<endl;
157  dfs::dfs_iterator itt2, endt2;
158  for (itt2 = search.begin(), endt2=search.end(); itt2 !=endt2; ++itt2)
159  {
160  cout<<*itt2<<" ";//<<"="<<search.dfs_num(*itt2)<<" ";
161  }
162  cout<<endl;
163  cout<<"Edge/Parton ordering result from DFS:"<<endl;
164  dfs::tree_edges_iterator itt, endt;
165  for (itt = search.tree_edges_begin(), endt=search.tree_edges_end(); itt !=endt; ++itt)
166  {
167  cout<<*itt;//<<endl;
168  }
169  cout<<endl;
170 
171  dfs::roots_iterator itt3, endt3;
172  cout<<"List of root nodes found in graph/shower : ";
173  for (itt3 = search.roots_begin(), endt3=search.roots_end(); itt3 !=endt3; ++itt3)
174  {
175  cout<<**itt3;
176  }
177  cout<<endl;
178  cout<<endl;
179 }
180 
181 
182 // -------------------------------------
183 
184 void Show()
185 {
186  ShowJetscapeBanner();
187  INFO_NICE;
188  INFO_NICE<<"------------------------------------";
189  INFO_NICE<<"| Reader Test JetScape Framework ... |";
190  INFO_NICE<<"------------------------------------";
191  INFO_NICE;
192 }
193 
194 //----------------------------------------------------------------------
196 
197 ostream & operator<<(ostream & ostr, const fjcore::PseudoJet & jet) {
198  if (jet == 0) {
199  ostr << " 0 ";
200  } else {
201  ostr << " pt = " << jet.pt()
202  << " m = " << jet.m()
203  << " y = " << jet.rap()
204  << " phi = " << jet.phi();
205  }
206  return ostr;
207 }
208 
209 
210 //----------------------------------------------------------------------