src/reader/JetScapeReader.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 "JetScapeReader.h"
0017 #include <sstream>
0018
0019 namespace Jetscape {
0020
0021 template <class T> JetScapeReader<T>::JetScapeReader():
0022   currentEvent{-1}
0023   , sigmaGen{-1}
0024   , sigmaErr{-1}
0025   , eventWeight{-1}
0026   , EventPlaneAngle{0.0}
0027 {
0028   VERBOSE(8);
0029 }
0030
0031 template <class T> JetScapeReader<T>::~JetScapeReader() { VERBOSE(8); }
0032
0033 template <class T> void JetScapeReader<T>::Clear() {
0034   nodeVec.clear();
0035   edgeVec.clear();
0036   //pShower->clear();//pShower=nullptr; //check ...
0037   pShowers.clear();
0038   hadrons.clear();
0039
0040   sigmaGen = -1;
0041   sigmaErr = -1;
0042   eventWeight = -1;
0043   EventPlaneAngle = 0.0;
0044 }
0045
0046 template <class T> void JetScapeReader<T>::AddNode(string s) {
0047   string token;
0048   //int counter=0;
0049   strT.set(s);
0050
0051   vector<string> vS;
0052
0053   while (!strT.done()) {
0054     token = strT.next();
0055     if (token.compare("V") != 0)
0056       vS.push_back(token);
0057   }
0058
0059   nodeVec.push_back(pShower->new_vertex(
0060       make_shared<Vertex>(stod(vS[1]), stod(vS[2]), stod(vS[3]), stod(vS[4]))));
0061 }
0062
0063 template <class T> void JetScapeReader<T>::AddEdge(string s) {
0064   if (nodeVec.size() > 1) {
0065     string token;
0066     //int counter=0;
0067     strT.set(s);
0068
0069     vector<string> vS;
0070
0071     while (!strT.done()) {
0072       token = strT.next();
0073       if (token.compare("P") != 0)
0074         vS.push_back(token);
0075     }
0076
0077     pShower->new_parton(
0078         nodeVec[stoi(vS[0])], nodeVec[stoi(vS[1])],
0079         make_shared<Parton>(
0080             stoi(vS[2]), stoi(vS[3]), stoi(vS[4]), stod(vS[5]), stod(vS[6]),
0081             stod(vS[7]),
0082             stod(
0083                 vS[8]))); // use different constructor wit true spatial posiiton ...
0084   } else
0085     JSWARN << "Node vector not filled, can not add edges/partons!";
0086 }
0087
0088 template <class T> void JetScapeReader<T>::AddHadron(string s) {
0089   string token;
0090   strT.set(s);
0091
0092   vector<string> vS;
0093   double x[4];
0094   x[0] = x[1] = x[2] = x[3] = 0.0;
0095   while (!strT.done()) {
0096     token = strT.next();
0097     if (token.compare("H") != 0)
0098       vS.push_back(token);
0099   }
0100   hadrons.push_back(make_shared<Hadron>(stoi(vS[1]), stoi(vS[2]), stoi(vS[3]),
0101                                         stod(vS[4]), stod(vS[5]), stod(vS[6]),
0102                                         stod(vS[7]), x));
0103 }
0104
0105 template <class T> void JetScapeReader<T>::Next() {
0106   if (currentEvent > 0)
0107     Clear();
0108
0109   //ReadEvent(currentPos);
0110   string line;
0111   string token;
0112
0113   JSINFO << "Current Event = " << currentEvent;
0114
0115   pShowers.push_back(make_shared<PartonShower>());
0116   pShower = pShowers[0];
0117   currentShower = 1;
0118
0119   int nodeZeroCounter = 0;
0120   std::string EPAngleStr = "EventPlaneAngle";
0121   while (getline(inFile, line)) {
0122     strT.set(line);
0123
0124     if (strT.isCommentEntry()) {
0125
0126       // Cross section
0127       if (line.find("sigmaGen") != std::string::npos) {
0128         std::stringstream data(line);
0129         std::string dummy;
0130         data >> dummy >> dummy >> sigmaGen;
0131         JSDEBUG << " sigma gen=" << sigmaGen;
0132       }
0133       // Cross section error
0134       if (line.find("sigmaErr") != std::string::npos) {
0135         std::stringstream data(line);
0136         std::string dummy;
0137         data >> dummy >> dummy >> sigmaErr;
0138         JSDEBUG << " sigma err=" << sigmaErr;
0139       }
0140       // Event weight
0141       if (line.find("weight") != std::string::npos) {
0142         std::stringstream data(line);
0143         std::string dummy;
0144         data >> dummy >> dummy >> eventWeight;
0145         JSDEBUG << " Event weight=" << eventWeight;
0146       }
0147       // EP angle
0148       if (line.find(EPAngleStr) != std::string::npos) {
0149         std::stringstream data(line);
0150         std::string dummy;
0151         data >> dummy >> dummy >> EventPlaneAngle;
0152         JSDEBUG << " EventPlaneAngle=" << EventPlaneAngle;
0153       }
0154       continue;
0155     }
0156
0157     if (strT.isEventEntry()) {
0158       int newEvent = stoi(strT.next());
0159       if (currentEvent != newEvent && currentEvent > -1) {
0160         currentEvent++;
0161         break;
0162       }
0163       currentEvent = newEvent;
0164       continue;
0165     }
0166
0167     // not an event header -- done?
0168     if (!strT.isGraphEntry())
0169       continue;
0170
0171     // node?
0172     if (strT.isNodeEntry()) {
0173       // catch starting node
0174       if (strT.isNodeZero()) {
0175         nodeZeroCounter++;
0176         if (nodeZeroCounter > currentShower) {
0177           nodeVec.clear();
0178           edgeVec.clear();
0179           pShowers.push_back(make_shared<PartonShower>());
0180           pShower = pShowers.back();
0181           currentShower++;
0182         }
0183       }
0184       AddNode(line);
0185       continue;
0186     }
0187
0188     // edge?
0189     if (strT.isEdgeEntry()) {
0190       AddEdge(line);
0191       continue;
0192     }
0193
0194     // rest is list entry == hadron entry
0195     // Some questionable nomenclature here - identifying all that begins with "[" as "GraphEntry"
0196     // oh well
0197     AddHadron(line);
0198   }
0199
0200   if (Finished())
0201     currentEvent++;
0202 }
0203
0204 template <class T>
0205 vector<fjcore::PseudoJet> JetScapeReader<T>::GetHadronsForFastJet() {
0206   vector<fjcore::PseudoJet> forFJ;
0207
0208   for (auto &h : hadrons) {
0209     forFJ.push_back(h->GetPseudoJet());
0210   }
0211
0212   return forFJ;
0213 }
0214
0215 template <class T> void JetScapeReader<T>::Init() {
0216   VERBOSE(8) << "Open Input File = " << file_name_in;
0217   JSINFO << "Open Input File = " << file_name_in;
0218
0219   inFile.open(file_name_in.c_str());
0220
0221   if (!inFile.good()) {
0222     JSWARN << "Corrupt input file!";
0223     exit(-1);
0224   } else
0225     JSINFO << "File opened";
0226
0227   currentEvent = 0;
0228 }
0229
0230 template class JetScapeReader<ifstream>;
0231
0232 #ifdef USE_GZIP
0233 template class JetScapeReader<igzstream>;
0234 #endif
0235
0236 } // end namespace Jetscape