JETSCAPE/examples/FinalStatePartons.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 depends on the header.
0062   std::string particleSeperator = " ";
0063   if (!writeHeaderV2) {
0064     particleSeperator = "\t";
0065   }
0066
0067   auto reader=make_shared<JetScapeReaderAscii>(argv[1]);
0068   std::ofstream dist_output (argv[2]); //Format is SN, PID, E, Px, Py, Pz, Eta, Phi
0069   int SN=0, TotalPartons=0;
0070   while (!reader->Finished())
0071   {
0072     reader->Next();
0073
0074     // cout<<"Analyze current event: "<<reader->GetCurrentEvent()<<endl;
0075     auto mShowers=reader->GetPartonShowers();
0076
0077     TotalPartons = 0;
0078     for (int i=0; i < mShowers.size(); i++)
0079     {
0080       TotalPartons = TotalPartons + mShowers[i]->GetFinalPartons().size();
0081     }
0082
0083     if(TotalPartons > 0)
0084     {
0085       ++SN;
0086       if (writeHeaderV2) {
0087         // NOTE: Needs consistent "\t" between all entries to simplify parsing later.
0088         dist_output << "#"
0089             << "\t" << "Event\t" << SN
0090             << "\t" << "weight\t" << reader->GetEventWeight()
0091             << "\t" << "EPangle\t" << reader->GetEventPlaneAngle()
0092             << "\t" << "N_partons\t" << TotalPartons
0093             << "\t" << "|"  // As a delimiter
0094             << "\t" << "N"
0095             << "\t" << "pid"
0096             << "\t" << "status"
0097             << "\t" << "E"
0098             << "\t" << "Px"
0099             << "\t" << "Py"
0100             << "\t" << "Pz"
0101             << "\n";
0102       }
0103       else {
0104         dist_output << "#" << "\t"
0105             << reader->GetEventPlaneAngle() << "\t"
0106             << "Event"
0107             << SN << "ID\t"
0108             << TotalPartons << "\t"
0109             << "pstat-EPx"  << "\t"
0110             << "Py"  << "\t"
0111             << "Pz"  << "\t"
0112             << "Eta" <<  "\t"<< "Phi" << "\n";
0113       }
0114
0115       for (int i=0; i < mShowers.size(); i++)
0116       {
0117         //cout<<" Analyze parton shower: "<<i<<endl;
0118         // Let's create a file
0119         for (int ipart = 0; ipart < mShowers[i]->GetFinalPartons().size(); ++ipart)
0120         {
0121           Parton p = *mShowers[i]->GetFinalPartons().at(ipart);
0122           //            if(abs(p.pid())!=5) continue;
0123
0124           dist_output << ipart
0125               << particleSeperator << p.pid()
0126               << particleSeperator << p.pstat()
0127               << particleSeperator << p.e()
0128               << particleSeperator << p.px()
0129               << particleSeperator << p.py()
0130               << particleSeperator << p.pz();
0131
0132           // v2 drops eta and phi, so only include it for v1
0133           if (!writeHeaderV2) {
0134               dist_output << particleSeperator << p.eta()
0135                   << particleSeperator << p.phi();
0136           }
0137
0138           // Finish up
0139           dist_output << "\n";
0140         }
0141       }
0142     }
0143   }
0144   // Write the final cross section and error if requested by using header v2
0145   if (writeHeaderV2) {
0146     // NOTE: Needs consistent "\t" between all entries to simplify parsing later.
0147     dist_output << "#"
0148         << "\t" << "sigmaGen\t" << reader->GetSigmaGen()
0149         << "\t" << "sigmaErr\t" << reader->GetSigmaErr()
0150         << "\n";
0151   }
0152   reader->Close();
0153 }