JETSCAPE/examples/readerTest.cc

0001 /*******************************************************************************
0002  * Copyright (c) The JETSCAPE Collaboration, 2018
0003  *
0004  * Modular, task-based framework for simulating all aspects of heavy-ion collisions
0005  *
0006  * For the list of contributors see AUTHORS.
0007  *
0008  * Report issues at https://github.com/JETSCAPE/JETSCAPE/issues
0009  *
0010  * or via email to bugs.jetscape@gmail.com
0011  *
0012  * Distributed under the GNU General Public License 3.0 (GPLv3 or later).
0013  * See COPYING for details.
0014  ******************************************************************************/
0015 // Reader test (focus on graph)
0016
0017 #include <iostream>
0018 #include <fstream>
0019 #include <memory>
0020 #include <chrono>
0021 #include <thread>
0022
0023 #include "gzstream.h"
0024 #include "PartonShower.h"
0025 #include "JetScapeLogger.h"
0026 #include "JetScapeReader.h"
0027 #include "JetScapeBanner.h"
0028 #include "fjcore.hh"
0029
0030 #include <GTL/dfs.h>
0031
0032 using namespace std;
0033 //using namespace fjcore;
0034
0035 using namespace Jetscape;
0036
0037 // -------------------------------------
0038
0039 // Forward declaration
0040 void Show();
0041 void AnalyzeGraph(shared_ptr<PartonShower> mS);
0042 ostream & operator<<(ostream & ostr, const fjcore::PseudoJet & jet);
0043
0044 // -------------------------------------
0045
0046 // Create a pdf of the shower graph:
0047 // Use with graphviz (on Mac: brew install graphviz --with-app)
0048 // in shell: dot GVfile.gv -Tpdf -o outputPDF.pdf
0049 // [or you can also use the GraphViz app for Mac Os X (in the "cellar" of homebrew)]
0050
0051 // -------------------------------------
0052
0053 int main(int argc, char** argv)
0054 {
0055   JetScapeLogger::Instance()->SetDebug(false);
0056   JetScapeLogger::Instance()->SetRemark(false);
0057   //SetVerboseLevel (9 a lot of additional debug output ...)
0058   //If you want to suppress it: use SetVerboseLevle(0) or max  SetVerboseLevle(9) or 10
0059   JetScapeLogger::Instance()->SetVerboseLevel(0);
0060
0061   cout<<endl;
0062   Show();
0063
0064   //Do some dummy jetfinding ...
0065   fjcore::JetDefinition jet_def(fjcore::antikt_algorithm, 0.7);
0066
0067   vector<shared_ptr<PartonShower>> mShowers;
0068
0069   //Directly with template: provide the relevant stream
0070   //auto reader=make_shared<JetScapeReader<ifstream>>("test_out.dat");
0071   //auto reader=make_shared<JetScapeReader<igzstream>>("test_out.dat.gz");
0072
0073   // Hide Template (see class declarations in reader/JetScapeReader.h) ...
0074   auto reader=make_shared<JetScapeReaderAscii>("test_out.dat");
0075   //auto reader=make_shared<JetScapeReaderAsciiGZ>("test_out.dat.gz");
0076
0077   // reads in multiple events and multiple shower per event
0078   // commentend out so that you get the dot graph file for the first shower in the first event
0079   // (add in and the file gets overriden)
0080   //while (!reader->Finished())
0081     {
0082       reader->Next();
0083
0084       cout<<"Analyze current event = "<<reader->GetCurrentEvent()<<endl;
0085       mShowers=reader->GetPartonShowers();
0086
0087       int finals = 0;
0088       for (int i=0;i<mShowers.size();i++)
0089     {
0090       cout<<" Analyze parton shower = "<<i<<endl;
0091
0092       mShowers[i]->PrintVertices();
0093       mShowers[i]->PrintPartons();
0094
0095       finals += mShowers[i]->GetFinalPartonsForFastJet().size();
0096
0097       fjcore::ClusterSequence cs(mShowers[i]->GetFinalPartonsForFastJet(), jet_def);
0098
0099       vector<fjcore::PseudoJet> jets = fjcore::sorted_by_pt(cs.inclusive_jets(2));
0100       cout<<endl;
0101       cout<<jet_def.description()<<endl;
0102       // Output of found jets ...
0103       //cout<<endl;
0104       for (int k=0;k<jets.size();k++)
0105         cout<<"Anti-kT jet "<<k<<" : "<<jets[k]<<endl;
0106       cout<<endl;
0107       cout<<"Shower initiating parton : "<<*(mShowers[i]->GetPartonAt(0))<<endl;
0108       cout<<endl;
0109
0110       AnalyzeGraph(mShowers[i]);
0111
0112       if (i==0)
0113         {
0114           mShowers[i]->SaveAsGV("my_test.gv");
0115           mShowers[i]->SaveAsGML("my_test.gml");
0116           mShowers[i]->SaveAsGraphML("my_test.graphml");
0117         }
0118
0119       // wait for 5s
0120       //std::this_thread::sleep_for(std::chrono::milliseconds(5000));
0121     }
0122       cout << " Found " << finals << " final state partons." << endl;
0123       auto hadrons = reader->GetHadrons();
0124       cout<<"Number of hadrons is: " << hadrons.size() << endl;
0125
0126       fjcore::ClusterSequence hcs(reader->GetHadronsForFastJet(), jet_def);
0127       vector<fjcore::PseudoJet> hjets = fjcore::sorted_by_pt(hcs.inclusive_jets(2));
0128       cout<<"AT HADRONIC LEVEL " << endl;
0129       for (int k=0;k<hjets.size();k++)
0130     cout<<"Anti-kT jet "<<k<<" : "<<hjets[k]<<endl;
0131
0132       // for(unsigned int i=0; i<hadrons.size(); i++) {
0133       //    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;
0134       // }
0135     }
0136
0137     reader->Close();
0138 }
0139
0140 // -------------------------------------
0141
0142 void AnalyzeGraph(shared_ptr<PartonShower> mS)
0143 {
0144   JSINFO<<"Some GTL graph/shower analysis/dfs search output:";
0145
0146   // quick and dirty ...
0147   dfs search;
0148   search.calc_comp_num(true);
0149   search.scan_whole_graph(true);
0150   search.start_node();// defaulted to first node ...
0151   search.run(*mS);
0152
0153   cout<<endl;
0154   cout<<"DFS graph search feature from GTL:"<<endl;
0155   cout<<"Number of Nodes reached from node 0 = "<< search.number_of_reached_nodes () <<endl;
0156   cout<<"Node/Vertex ordering result from DFS:"<<endl;
0157   dfs::dfs_iterator itt2, endt2;
0158   for (itt2 = search.begin(), endt2=search.end(); itt2 !=endt2; ++itt2)
0159     {
0160       cout<<*itt2<<" ";//<<"="<<search.dfs_num(*itt2)<<" ";
0161     }
0162   cout<<endl;
0163   cout<<"Edge/Parton ordering result from DFS:"<<endl;
0164   dfs::tree_edges_iterator itt, endt;
0165   for (itt = search.tree_edges_begin(), endt=search.tree_edges_end(); itt !=endt; ++itt)
0166     {
0167         cout<<*itt;//<<endl;
0168     }
0169   cout<<endl;
0170
0171   dfs::roots_iterator itt3, endt3;
0172   cout<<"List of root nodes found in graph/shower : ";
0173   for (itt3 = search.roots_begin(), endt3=search.roots_end(); itt3 !=endt3; ++itt3)
0174     {
0175       cout<<**itt3;
0176     }
0177   cout<<endl;
0178   cout<<endl;
0179 }
0180
0181
0182 // -------------------------------------
0183
0184 void Show()
0185 {
0186   ShowJetscapeBanner();
0187   INFO_NICE;
0188   INFO_NICE<<"------------------------------------";
0189   INFO_NICE<<"| Reader Test JetScape Framework ... |";
0190   INFO_NICE<<"------------------------------------";
0191   INFO_NICE;
0192 }
0193
0194 //----------------------------------------------------------------------
0195 /// overloaded jet info output
0196
0197 ostream & operator<<(ostream & ostr, const fjcore::PseudoJet & jet) {
0198   if (jet == 0) {
0199     ostr << " 0 ";
0200   } else {
0201     ostr << " pt = " << jet.pt()
0202          << " m = " << jet.m()
0203          << " y = " << jet.rap()
0204          << " phi = " << jet.phi();
0205   }
0206   return ostr;
0207 }
0208
0209
0210 //----------------------------------------------------------------------