JETSCAPE/examples/FinalStateHadrons.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
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 using namespace std;
0032 using namespace fjcore;
0033
0034 using namespace Jetscape;
0035
0036 // You could overload here and then simply use ofstream << p;
0037 // ostream & operator<<(ostream & ostr, const fjcore::PseudoJet & jet);
0038
0039
0040 // -------------------------------------
0041
0042 int main(int argc, char** argv)
0043 {
0044
0045   // JetScapeLogger::Instance()->SetInfo(false);
0046   JetScapeLogger::Instance()->SetDebug(false);
0047   JetScapeLogger::Instance()->SetRemark(false);
0048   // //SetVerboseLevel (9 a lot of additional debug output ...)
0049   // //If you want to suppress it: use SetVerboseLevle(0) or max  SetVerboseLevel(9) or 10
0050   JetScapeLogger::Instance()->SetVerboseLevel(0);
0051
0052   // Whether to write the new header (ie. v2), including xsec info.
0053   // To enable, pass anything as the third argument to enable this option.
0054   // Default: disabled.
0055   bool writeHeaderV2 = false;
0056   if (argc > 3) {
0057     writeHeaderV2 = static_cast<bool>(argv[3]);
0058     std::cout << "NOTE: Writing header v2, and final cross section and error at EOF.\n";
0059   }
0060
0061   // The seperator between particles _does not_ depend on the header for final state hadrons
0062   std::string particleSeperator = " ";
0063
0064   auto reader = make_shared<JetScapeReaderAscii>(argv[1]);
0065   std::ofstream dist_output (argv[2]); //Format is SN, PID, E, Px, Py, Pz, Eta, Phi
0066   vector<shared_ptr<Hadron>> hadrons;
0067   int SN=0;
0068   while (!reader->Finished())
0069   {
0070     reader->Next();
0071
0072     // cout<<"Analyze current event: "<<reader->GetCurrentEvent()<<endl;
0073
0074     //dist_output<<"Event "<< reader->GetCurrentEvent()+1<<endl;
0075     hadrons = reader->GetHadrons();
0076     cout<<"Number of hadrons is: " << hadrons.size() << endl;
0077
0078     if (hadrons.size() > 0)
0079     {
0080       ++SN;
0081       if (writeHeaderV2) {
0082         // NOTE: Needs consistent "\t" between all entries to simplify parsing later.
0083         dist_output << "#"
0084             << "\t" << "Event\t" << SN
0085             << "\t" << "weight\t" << reader->GetEventWeight()
0086             << "\t" << "EPangle\t" << reader->GetEventPlaneAngle()
0087             << "\t" << "N_hadrons\t" << hadrons.size()
0088             << "\t" << "|"  // As a delimiter
0089             << "\t" << "N"
0090             << "\t" << "pid"
0091             << "\t" << "status"
0092             << "\t" << "E"
0093             << "\t" << "Px"
0094             << "\t" << "Py"
0095             << "\t" << "Pz"
0096             << "\n";
0097       }
0098       else {
0099         dist_output << "#" << "\t"
0100             << reader->GetEventPlaneAngle() << "\t"
0101             << "Event"
0102             << SN << "ID\t"
0103             << hadrons.size() << "\t"
0104             << "pstat-EPx"   << "\t"
0105             << "Py"  << "\t"
0106             << "Pz"  << "\t"
0107             << "Eta" <<  "\t"<< "Phi" << "\n";
0108       }
0109
0110       for (unsigned int i=0; i<hadrons.size(); i++)
0111       {
0112         dist_output << i
0113             << particleSeperator << hadrons[i].get()->pid()
0114             << particleSeperator << hadrons[i].get()->pstat()
0115             << particleSeperator << hadrons[i].get()->e()
0116             << particleSeperator << hadrons[i].get()->px()
0117             << particleSeperator << hadrons[i].get()->py()
0118             << particleSeperator << hadrons[i].get()->pz();
0119
0120         // v2 drops eta and phi, so only include it for v1
0121         if (!writeHeaderV2) {
0122             dist_output << particleSeperator << hadrons[i].get()->eta()
0123                 << particleSeperator << hadrons[i].get()->phi();
0124         }
0125
0126         // Finish up
0127         dist_output << "\n";
0128       }
0129     }
0130   }
0131   // Write the final cross section and error if requested by using header v2
0132   if (writeHeaderV2) {
0133     // NOTE: Needs consistent "\t" between all entries to simplify parsing later.
0134     dist_output << "#"
0135         << "\t" << "sigmaGen\t" << reader->GetSigmaGen()
0136         << "\t" << "sigmaErr\t" << reader->GetSigmaErr()
0137         << "\n";
0138   }
0139   reader->Close();
0140 }