JETSCAPE/jail/pwg2_reader.cxx

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
0016 #include <iostream>
0017 #include <fstream>
0018 #include <memory>
0019 #include <chrono>
0020 #include <thread>
0021
0022 #include "gzstream.h"
0023 #include "PartonShower.h"
0024 #include "JetScapeLogger.h"
0025 #include "JetScapeReader.h"
0026 #include "JetScapeBanner.h"
0027 #include "fjcore.hh"
0028
0029 #include <GTL/dfs.h>
0030
0031 #include <TH1.h>
0032 #include <TFile.h>
0033
0034 using namespace std;
0035 using namespace Jetscape;
0036
0037 using namespace fjcore;
0038
0039 // -------------------------------------
0040
0041 // Forward declaration
0042 void AnalyzeGraph(shared_ptr<PartonShower> mS);
0043 ostream & operator<<(ostream & ostr, const PseudoJet & jet);
0044
0045 // -------------------------------------
0046
0047 // Create a pdf of the shower graph:
0048 // Use with graphviz (on Mac: brew install graphviz --with-app)
0049 // in shell: dot GVfile.gv -Tpdf -o outputPDF.pdf
0050 // [or you can also use the GraphViz app for Mac Os X (in the "cellar" of homebrew)]
0051
0052 // -------------------------------------
0053
0054 int main(int argc, char** argv)
0055 {
0056
0057   // Parse arguments
0058   string inname  = "test_out.dat.gz";
0059   string outname = "test_out.root";
0060   string filetype = ".dat.gz";
0061   std::size_t pos=0;
0062   switch ( argc ){
0063     break;
0064   case 3:
0065     inname = argv[1];
0066     outname = argv[2];
0067     break;
0068   case 2:
0069     inname = argv[1];
0070     outname = inname;
0071     pos = outname.find( filetype);
0072     if ( pos==std::string::npos ){
0073       throw std::runtime_error ("Please use a filename ending in .dat.gz for an automatically generated output name");
0074       return -1;
0075     }
0076     outname.replace(pos, filetype.length(), ".root");
0077     break;
0078   case 1:
0079   case 0:
0080     break;
0081   default:
0082     cout << "Usage: readerTestWithRoot [inputname] [outputname]" << endl;
0083     return -1;
0084   }
0085
0086   JetScapeLogger::Instance()->SetInfo(false);
0087   JetScapeLogger::Instance()->SetDebug(false);
0088   JetScapeLogger::Instance()->SetRemark(false);
0089   JetScapeLogger::Instance()->SetVerboseLevel(0);
0090
0091   cout<<endl;
0092
0093   // VERY crude way to pull event information that
0094   // resides at the foot of the file in a comment section
0095   // (for technical reasons);
0096   igzstream inFile;
0097   inFile.open(inname.c_str());
0098   string line;
0099   int nAccepted=-1;
0100   double sigmaGen=-1;
0101   string cs="# ";
0102   string tofind;
0103   while (getline(inFile,line)) {
0104     if ( line.compare(0, cs.length(), cs) != 0) continue;
0105
0106     tofind = "# nAccepted = ";
0107     pos = line.find( tofind );
0108     if ( pos!=std::string::npos )      nAccepted=atoi( string ( line.begin()+pos+tofind.length(), line.end()).c_str());
0109
0110     tofind = "# sigmaGen  = ";
0111     pos = line.find( tofind );
0112     if ( pos!=std::string::npos )      sigmaGen=atof( string ( line.begin()+pos+tofind.length(), line.end()).c_str());
0113   }
0114   if ( sigmaGen<0 || nAccepted<0 ){
0115     throw std::runtime_error("Couldn't find sigmaGen and nAccepted in the footer.");
0116     return -1;
0117   }
0118
0119   INFO << "Using sigmaGen = " << sigmaGen << " and nAccepted = " << nAccepted;
0120
0121   double weight = sigmaGen / nAccepted;
0122   // weight =1;
0123
0124   // Simple jetfinding ...
0125   float R=0.4;
0126   JetDefinition jet_def(antikt_algorithm, R);
0127
0128   // Constituent selectors
0129   // ---------------------
0130   float max_track_rap = 2;
0131   Selector select_track_rap = SelectorAbsRapMax(max_track_rap);
0132   Selector select_pt        = SelectorPtMin( 0.2 );
0133   Selector select_cons      = select_track_rap * select_pt;
0134
0135   // jet selectors
0136   Selector select_jet_eta = SelectorAbsEtaMax( max_track_rap-R );
0137   Selector select_jet_pt  = SelectorPtRange( 10,1000 );
0138   Selector select_jet     = select_jet_eta * select_jet_pt;
0139
0140
0141   // dijet cuts
0142   float ptlead = 40;
0143   float ptsublead = 20;
0144
0145   auto reader=make_shared<JetScapeReaderAsciiGZ>(inname);
0146
0147   // Open output root file
0148   TFile* fout = new TFile( outname.data(), "RECREATE");
0149
0150   // turn on error bars by default
0151   TH1::SetDefaultSumw2();
0152
0153   // Histograms
0154   TH1D* partonpt = new TH1D ("partonpt","Final Parton p_{T};p_{T} [GeV/c]",120, 0, 120);
0155   TH1D* partoneta = new TH1D ("partoneta","Final Parton #eta; #eta",100, -5, 5);
0156   TH1D* hardpartonpt = new TH1D ("hardpartonpt","Hard Parton p_{T};p_{T} [GeV/c]",80, 0, 40);
0157   TH1D* hardpartoneta = new TH1D ("hardpartoneta","Hard Parton #eta; #eta",100, -5, 5);
0158
0159   TH1D* leadjetpt = new TH1D ("leadjetpt","Leading Jet p_{T};p_{T}^{jet} [GeV/c]",120, 0, 120);
0160   TH1D* subleadjetpt = new TH1D ("subleadjetpt","Sub-Leading Jet p_{T};p_{T}^{jet} [GeV/c]",120, 0, 120);
0161
0162   TH1D* aj = new TH1D ("aj","A_{J};", 40, -1, 1 );
0163
0164   TH1D* jetpt = new TH1D ("jetpt","Jet p_{T};p_{T}^{jet} [GeV/c]",120,0,120);
0165   TH1D* jeteta = new TH1D ("jeteta","Jet #eta;#eta^{jet} [GeV/c]",100, -3, 3);
0166
0167   // Histogram for jet fragmentation function
0168   TH1D* jetz = new TH1D ("jetz","Jet z;p_{T}^{constituent} / p_{T}^{jet}",100, 0, 1);
0169
0170   TH1D* dR = new TH1D ("dR","#DeltaR;counts",50, 0, 1);
0171
0172   // reads in multiple events and multiple shower per event
0173   vector<shared_ptr<PartonShower>> mShowers;
0174
0175   long event = 0;
0176   while (!reader->Finished()){
0177     reader->Next();
0178     if ( !(event%1000) ) cout << "Working on " << event << endl;
0179     event++;
0180
0181     // INFO<<"Analyze current event = "<<reader->GetCurrentEvent();
0182     mShowers=reader->GetPartonShowers();
0183
0184     // Collect jets from all showers
0185     vector<PseudoJet> jets;
0186     for (int i=0;i<mShowers.size();i++){
0187       // Actual jetfinding
0188       vector<PseudoJet> particles = select_cons(mShowers[i]->GetFinalPartonsForFastJet());
0189       ClusterSequence cs(particles, jet_def);
0190
0191       // mShowers[i]->PrintVertices();
0192       // mShowers[i]->PrintPartons();
0193
0194       // inital hard partons
0195       hardpartonpt->Fill ( mShowers[i]->GetPartonAt(0)->pt(), weight );
0196       hardpartoneta->Fill ( mShowers[i]->GetPartonAt(0)->eta(), weight );
0197
0198       // final particles
0199       for ( auto &p : mShowers[i]->GetFinalPartonsForFastJet() ){
0200     partonpt->Fill ( p.pt(), weight );
0201     partoneta->Fill ( p.eta(), weight );
0202       }
0203
0204       // jets
0205       for (auto j : select_jet( cs.inclusive_jets()) ){
0206     // additional cuts?
0207     jets.push_back(j);
0208
0209     // save some constituent information
0210     for ( auto &consts : j.constituents() ){
0211       // if ( consts.pt() > 40 ) {
0212       //   cout << "---" << endl;
0213       //   cout << consts << endl;
0214       //   // cout << " --> " << *(mShowers[i]->GetPartonAt(0)) << endl;
0215       //   cout << " --> shower initiator at pt = " << mShowers[i]->GetPartonAt(0)->pt()
0216       //       << " eta = " << mShowers[i]->GetPartonAt(0)->eta()  << endl;
0217       //   // for ( auto &parts : mShowers[i]->GetFinalPartons() ){
0218       //   //   cout << " ----> " << *parts << endl;
0219       //   //   cout << " ----> " << parts->pt() << endl;
0220       //   // }
0221
0222       // }
0223       jetz->Fill( consts.pt() / j.pt(), weight );
0224       dR->Fill( consts.delta_R(j), weight );
0225     }
0226
0227       }
0228
0229
0230     }
0231
0232     // analyze jets
0233     if ( jets.size()==0 ) continue;
0234
0235     jets = sorted_by_pt( jets );
0236
0237     // misc information
0238     for (auto j : jets){
0239       // Fill histograms
0240       jetpt->Fill( j.pt(), weight );
0241       jeteta->Fill( j.eta(), weight );
0242
0243     }
0244
0245     // Dijets
0246     if ( jets.size()<2 ) continue;
0247     PseudoJet lead = jets.at(0);
0248     PseudoJet sublead = jets.at(1);
0249
0250     if ( lead.pt()<ptlead || sublead.pt()<ptsublead ) continue;
0251
0252     // back to back?
0253     float dphi = lead.delta_phi_to( sublead ); // always in -pi,pi
0254     if ( dphi < -Jetscape::pi || dphi > Jetscape::pi ) {
0255       cerr << " dphi = " << dphi << endl;
0256       return -1;
0257     }
0258
0259     if ( dphi > -Jetscape::pi+R && dphi < Jetscape::pi-R ) continue;
0260
0261     // Fill histos
0262     leadjetpt->Fill( lead.pt(), weight );
0263     subleadjetpt->Fill( sublead.pt(), weight );
0264     aj->Fill( (lead.pt()-sublead.pt()) / (lead.pt()+sublead.pt()), weight);
0265
0266
0267   }
0268   reader->Close();
0269   fout->Write();
0270
0271   INFO_NICE<<"Finished!";
0272   INFO << "Read from " <<  inname;
0273   INFO << "Wrote to " <<  outname;
0274   cout << endl;
0275   return 0;
0276
0277 }
0278
0279 // -------------------------------------
0280
0281 void AnalyzeGraph(shared_ptr<PartonShower> mS)
0282 {
0283   INFO<<"Some GTL graph/shower analysis/dfs search output:";
0284
0285   dfs search;
0286   search.calc_comp_num(true);
0287   search.scan_whole_graph(true);
0288   search.start_node();// defaulted to first node ...
0289   search.run(*mS);
0290
0291   cout<<endl;
0292   cout<<"DFS graph search feature from GTL:"<<endl;
0293   cout<<"Number of Nodes reached from node 0 = "<< search.number_of_reached_nodes () <<endl;
0294   cout<<"Node/Vertex ordering result from DFS:"<<endl;
0295   dfs::dfs_iterator itt2, endt2;
0296   for (itt2 = search.begin(), endt2=search.end(); itt2 !=endt2; ++itt2)
0297     {
0298       cout<<*itt2<<" ";//<<"="<<search.dfs_num(*itt2)<<" ";
0299     }
0300   cout<<endl;
0301
0302   cout<<"Edge/Parton ordering result from DFS:"<<endl;
0303   dfs::tree_edges_iterator itt, endt;
0304   for (itt = search.tree_edges_begin(), endt=search.tree_edges_end(); itt !=endt; ++itt)
0305     {
0306       cout<<*itt;//<<endl;
0307     }
0308   cout<<endl;
0309
0310   dfs::roots_iterator itt3, endt3;
0311   cout<<"List of root nodes found in graph/shower : ";
0312   for (itt3 = search.roots_begin(), endt3=search.roots_end(); itt3 !=endt3; ++itt3)
0313     {
0314       cout<<**itt3;
0315     }
0316   cout<<endl;
0317   cout<<endl;
0318 }
0319
0320 //----------------------------------------------------------------------
0321 /// overloaded jet info output
0322 ostream & operator<<(ostream & ostr, const PseudoJet & jet) {
0323   if (jet == 0) {
0324     ostr << " 0 ";
0325   } else {
0326     ostr << " pt = " << jet.pt()
0327          << " m = " << jet.m()
0328          << " y = " << jet.rap()
0329          << " phi = " << jet.phi();
0330   }
0331   return ostr;
0332 }
0333
0334
0335 //----------------------------------------------------------------------