src/framework/PartonShower.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 "PartonShower.h"
0017 #include <sstream>
0018 #include <fstream>
0019 #include <iomanip>
0020 #include "MakeUniqueHelper.h"
0021
0022 using std::setprecision;
0023 using std::fixed;
0024 using std::to_string;
0025 using std::ostringstream;
0026 using std::stringstream;
0027
0028 namespace Jetscape {
0029
0030 PartonShower::PartonShower() : graph() { VERBOSESHOWER(8); }
0031
0032 node PartonShower::new_vertex(shared_ptr<Vertex> v) {
0033   node n = graph::new_node();
0034   vMap[n] = v;
0035   return n;
0036 }
0037
0038 int PartonShower::new_parton(node s, node t, shared_ptr<Parton> p) {
0039   edge e = graph::new_edge(s, t);
0040   pMap[e] = p;
0041   return e.id();
0042 }
0043
0044 /*
0045 void PartonShower::FillPartonVec()
0046 {
0047   VERBOSE(8);
0048
0049   edge_iterator eIt, eEnd;
0050
0051   for (eIt = edges_begin(), eEnd = edges_end(); eIt != eEnd; ++eIt)
0052     pVec.push_back(pMap[*eIt]);
0053 }
0054
0055 void PartonShower::FillVertexVec()
0056 {
0057   VERBOSE(8);
0058
0059   node_iterator nIt, nEnd;
0060
0061   for (nIt = nodes_begin(), nEnd = nodes_end(); nIt != nEnd; ++nIt)
0062     vVec.push_back(vMap[*nIt]);
0063 }
0064 */
0065
0066 /*
0067 shared_ptr<Parton> PartonShower::GetPartonAt(int i)
0068 {
0069   if (pVec.size()==0)
0070     FillPartonVec();
0071
0072   return pVec[i].lock();
0073 }
0074
0075 shared_ptr<Vertex> PartonShower::GetVertexAt(int i)
0076 {
0077   if (vVec.size()==0)
0078     FillVertexVec();
0079
0080   return vVec[i].lock();
0081 }
0082 */
0083
0084 /*
0085 unique_ptr<Parton> PartonShower::GetPartonAt(int i)
0086 {
0087   if (pVec.size()==0)
0088     FillPartonVec();
0089
0090   return make_unique<Parton>(*(pVec[i].lock()));
0091 }
0092
0093 unique_ptr<Vertex> PartonShower::GetVertexAt(int i)
0094 {
0095   if (vVec.size()==0)
0096     FillVertexVec();
0097
0098   return make_unique<Vertex>(*(vVec[i].lock()));
0099 }
0100 */
0101 /*
0102 void PartonShower::CreateMaps()
0103 {
0104   VERBOSESHOWER(8);
0105   edge_iterator eIt, eEnd;
0106   for (eIt = edges_begin(), eEnd = edges_end(); eIt != eEnd; ++eIt)
0107     {
0108       pToEdgeMap[pMap[*eIt]]=*eIt;
0109     }
0110
0111   node_iterator nIt, nEnd;
0112   for (nIt = nodes_begin(), nEnd = nodes_end(); nIt != nEnd; ++nIt)
0113     {
0114       //vToNodeMap[vMap[*nIt]]=*nIt;
0115     }
0116 }
0117 */
0118
0119 vector<shared_ptr<Parton>> PartonShower::GetFinalPartons() {
0120   //VERBOSESHOWER(8)<<pFinal.size();
0121   if (pFinal.size() == 0) {
0122     edge_iterator eIt, eEnd;
0123     for (eIt = edges_begin(), eEnd = edges_end(); eIt != eEnd; ++eIt) {
0124       if (eIt->target().outdeg() < 1) {
0125         if (pMap[*eIt]->pstat() > -10) {
0126           pFinal.push_back(pMap[*eIt]);
0127         }
0128         // DEBUG
0129         //cout<<eIt->target()<<endl;
0130       }
0131     }
0132     return pFinal;
0133   } else {
0134     return pFinal;
0135   }
0136 }
0137
0138 vector<fjcore::PseudoJet> PartonShower::GetFinalPartonsForFastJet() {
0139   vector<shared_ptr<Parton>> mP =
0140       GetFinalPartons(); // to ensure that pFinal is filled
0141
0142   vector<fjcore::PseudoJet> forFJ;
0143
0144   for (int i = 0; i < mP.size(); i++)
0145     forFJ.push_back((*mP[i]).GetPseudoJet());
0146
0147   return forFJ;
0148 }
0149
0150 int PartonShower::GetNumberOfParents(int n) {
0151   return GetEdgeAt(n).source().indeg();
0152 }
0153
0154 int PartonShower::GetNumberOfChilds(int n) {
0155   return GetEdgeAt(n).target().outdeg();
0156 }
0157
0158 edge PartonShower::GetEdgeAt(int n) {
0159   edge_iterator eIt;
0160   eIt = edges_begin();
0161   advance(eIt, n);
0162   return *eIt;
0163 }
0164
0165 node PartonShower::GetNodeAt(int n) {
0166   node_iterator nIt;
0167   nIt = nodes_begin();
0168   advance(nIt, n);
0169   return *nIt;
0170 }
0171
0172 shared_ptr<Parton> PartonShower::GetPartonAt(int n) {
0173   edge_iterator eIt;
0174   eIt = edges_begin();
0175   advance(eIt, n);
0176   return GetParton(*eIt);
0177 }
0178
0179 shared_ptr<Vertex> PartonShower::GetVertexAt(int n) {
0180   node_iterator nIt;
0181   nIt = nodes_begin();
0182   advance(nIt, n);
0183   return GetVertex(*nIt);
0184 }
0185
0186 PartonShower::~PartonShower() {
0187   VERBOSESHOWER(8);
0188   pFinal.clear(); //pVec.clear();vVec.clear();
0189 }
0190
0191 void PartonShower::save_node_info_handler(ostream *o, node n) const {
0192   *o << "label "
0193      << "\"" << n.id() << "(" << fixed << setprecision(2) << vMap[n]->x_in().t()
0194      << ")\"" << endl;
0195   *o << "x " << vMap[n]->x_in().x() << endl;
0196   *o << "y " << vMap[n]->x_in().y() << endl;
0197   *o << "z " << vMap[n]->x_in().z() << endl;
0198   *o << "t " << vMap[n]->x_in().t() << endl;
0199 }
0200
0201 void PartonShower::save_edge_info_handler(ostream *o, edge e) const {
0202   *o << "label "
0203      << "\"(" << fixed << setprecision(2) << pMap[e]->pt() << ")\"" << endl;
0204   *o << "plabel " << pMap[e]->plabel() << endl;
0205   *o << "pid " << pMap[e]->pid() << endl;
0206   *o << "pstat " << pMap[e]->pstat() << endl;
0207   *o << "pT " << pMap[e]->pt() << endl;
0208   *o << "eta " << pMap[e]->eta() << endl;
0209   *o << "phi " << pMap[e]->phi() << endl;
0210   *o << "E " << pMap[e]->e() << endl;
0211 }
0212
0213 void PartonShower::pre_clear_handler() {
0214   VERBOSESHOWER(8);
0215   edge_iterator eIt, eEnd;
0216   for (eIt = edges_begin(), eEnd = edges_end(); eIt != eEnd; ++eIt) {
0217     pMap[*eIt] = nullptr;
0218   }
0219
0220   node_iterator nIt, nEnd;
0221   for (nIt = nodes_begin(), nEnd = nodes_end(); nIt != nEnd; ++nIt) {
0222     vMap[*nIt] = nullptr;
0223   }
0224 }
0225
0226 void PartonShower::PrintNodes(bool verbose) {
0227   node_iterator nIt, nEnd;
0228   ostringstream os;
0229
0230   if (verbose && JetScapeLogger::Instance()->GetVerboseLevel() > 8) {
0231     for (nIt = nodes_begin(), nEnd = nodes_end(); nIt != nEnd; ++nIt)
0232       os << *nIt << "=" << vMap[*nIt]->x_in().t() << " ";
0233     VERBOSESHOWER(8) << os.str();
0234   }
0235
0236   if (!verbose) {
0237     for (nIt = nodes_begin(), nEnd = nodes_end(); nIt != nEnd; ++nIt)
0238       os << *nIt << "=" << vMap[*nIt]->x_in().t() << " ";
0239     cout << "Vertex list : " << os.str() << endl;
0240   }
0241 }
0242
0243 void PartonShower::PrintEdges(bool verbose) {
0244   edge_iterator eIt, eEnd;
0245   ostringstream os;
0246
0247   if (verbose && JetScapeLogger::Instance()->GetVerboseLevel() > 8) {
0248     for (eIt = edges_begin(), eEnd = edges_end(); eIt != eEnd; ++eIt)
0249       os << *eIt << "=" << pMap[*eIt]->pt() << " ";
0250     VERBOSESHOWER(8) << os.str();
0251   }
0252
0253   if (!verbose) {
0254     for (eIt = edges_begin(), eEnd = edges_end(); eIt != eEnd; ++eIt)
0255       os << *eIt << "=" << pMap[*eIt]->pt() << " ";
0256     cout << "Parton list : " << os.str() << endl;
0257   }
0258 }
0259
0260 // To be extended to store all infos like this in GML
0261 // or store vectore of partons and vertecies with relevant graph infos
0262 // to construct the graph later .... (TBD)
0263 // These handlers are needed if one wants to use load and GML files ...
0264 // Obsolete with new JetScape reader ...
0265
0266 void PartonShower::load_edge_info_handler(edge e, GML_pair *read) {
0267   VERBOSESHOWER(8) << "Load edge ... " << e;
0268   struct GML_pair *tmp = read;
0269
0270   double pT, eta, phi, E;
0271   int plabel, pid, pstat;
0272
0273   while (tmp) {
0274     //printf ("*KEY* : %s \n", tmp->key);
0275     if (((string)(tmp->key)).find("pT") < 1)
0276       pT = tmp->value.floating;
0277     if (((string)(tmp->key)).find("eta") < 1)
0278       eta = tmp->value.floating;
0279     if (((string)(tmp->key)).find("phi") < 1)
0280       phi = tmp->value.floating;
0281     if (((string)(tmp->key)).find("E") < 1)
0282       E = tmp->value.floating;
0283     if (((string)(tmp->key)).find("plabel") < 1)
0284       plabel = tmp->value.integer;
0285     if (((string)(tmp->key)).find("pid") < 1)
0286       pid = tmp->value.integer;
0287     if (((string)(tmp->key)).find("pstat") < 1)
0288       pstat = tmp->value.integer;
0289
0290     tmp = tmp->next;
0291   }
0292
0293   pMap[e] = make_shared<Parton>(plabel, pid, pstat, pT, eta, phi, E);
0294 }
0295
0296 void PartonShower::load_node_info_handler(node n, GML_pair *read) {
0297   VERBOSESHOWER(8) << "Load node ... " << n;
0298   struct GML_pair *tmp = read;
0299
0300   double x = 0;
0301   double y = 0;
0302   double z = 0;
0303   double t = 0;
0304
0305   while (tmp) {
0306     //printf ("*KEY* : %s %f \n", tmp->key, tmp->value.floating);
0307     if (((string)(tmp->key)).find("x") < 1)
0308       x = tmp->value.floating;
0309     if (((string)(tmp->key)).find("y") < 1)
0310       y = tmp->value.floating;
0311     if (((string)(tmp->key)).find("z") < 1)
0312       z = tmp->value.floating;
0313     if (((string)(tmp->key)).find("t") < 1)
0314       t = tmp->value.floating;
0315
0316     tmp = tmp->next;
0317   }
0318
0319   vMap[n] = make_shared<Vertex>(x, y, z, t);
0320 }
0321
0322 // use with graphviz (on Mac: brew install graphviz --with-app)
0323 // dot GVfile.gv -Tpdf -o outputPDF.pdf
0324
0325 void PartonShower::SaveAsGV(string fName) {
0326   ofstream gv;
0327   gv.open(fName.c_str());
0328
0329   // Simple directed graph left->right in dot/gv format for usage with graphviz ...
0330   // nodes show (time) and arrows (pT)
0331   gv << "digraph \"graph\" {" << endl;
0332   gv << endl;
0333   gv << "rankdir=\"LR\";" << endl;
0334   gv << "node [shape=plaintext, fontsize=11];"
0335      << endl; //, shape=circle]; //plaintext];
0336   gv << "edge [fontsize=10];" << endl;
0337   gv << endl;
0338   //gv<<"0 -> 1"<<endl;
0339   node_iterator nIt, nEnd;
0340
0341   int n = 0;
0342   string label;
0343   string label2;
0344
0345   for (nIt = nodes_begin(), nEnd = nodes_end(); nIt != nEnd; ++nIt) {
0346     label = ("[label=\"" + to_string(n) + "(");
0347
0348     label2 = ")\"];";
0349     stringstream stream;
0350
0351     stream << fixed << setprecision(2) << (vMap[*nIt]->x_in().t());
0352     gv << n << " " << label << stream.str() << label2 << endl;
0353
0354     n++;
0355   }
0356
0357   gv << endl;
0358
0359   edge_iterator eIt, eEnd;
0360
0361   for (eIt = edges_begin(), eEnd = edges_end(); eIt != eEnd; ++eIt) {
0362     //label = ("[label=\"(");
0363     if ((pMap[*eIt]->pstat()) == -13) {
0364       // missing energy-momentum
0365       label = ("[style=\"dotted\" color=\"red\"label=\"(");
0366     } else if ((pMap[*eIt]->pstat()) == -11) {
0367       // liquefied partons
0368       label = ("[color=\"red\"label=\"(");
0369     } else if ((pMap[*eIt]->pstat()) == -17) {
0370       // thermal partons draw from the medium (negative partons)
0371       label = ("[color=\"darkorange\"label=\"(");
0372     } else if ((pMap[*eIt]->pstat()) == -1) {
0373       // thermal partons draw from the medium (negative partons)
0374       label = ("[color=\"green\"label=\"(");
0375     } else if ((pMap[*eIt]->t()) < 4.0) {
0376       // small virtuality parton
0377       label = ("[color=\"blue\"label=\"(");
0378     } else {
0379       label = ("[label=\"(");
0380     }
0381     label2 = ")\"];";
0382     stringstream stream;
0383     if ((pMap[*eIt]->pstat()) == -13) {
0384       stream << std::scientific << setprecision(1) << (pMap[*eIt]->e()) << ","
0385              << (pMap[*eIt]->pt()) << "," << (pMap[*eIt]->t()) << ","
0386              << (pMap[*eIt]->pid());
0387     } else {
0388       stream << fixed << setprecision(2) << (pMap[*eIt]->e()) << ","
0389              << (pMap[*eIt]->pt()) << "," << (pMap[*eIt]->t()) << ","
0390              << (pMap[*eIt]->pid());
0391     }
0392
0393     gv << (to_string(eIt->source().id()) + "->" + to_string(eIt->target().id()))
0394        << " " << label << stream.str() << label2 << endl;
0395   }
0396   gv << endl;
0397   gv << "}" << endl;
0398   gv.close();
0399 }
0400
0401 void PartonShower::SaveAsGraphML(string fName) {
0402   // Think about using tinyxml2 in future (if needed ...)
0403
0404   ofstream g;
0405   g.open(fName.c_str());
0406
0407   g << "<?xml version=\"1.0\" encoding=\"UTF-8\"?>" << endl;
0408   g << "<graphml xmlns=\"http://graphml.graphdrawing.org/xmlns\"" << endl;
0409   g << "xmlns:xsi=\"http://www.w3.org/2001/XMLSchema-instance\"" << endl;
0410   g << "xsi:schemaLocation=\"http://graphml.graphdrawing.org/xmlns" << endl;
0411   g << " http://graphml.graphdrawing.org/xmlns/1.0/graphml.xsd\">" << endl;
0412
0413   g << "<key id=\"nlabel\" for=\"node\" attr.name=\"label\" "
0414        "attr.type=\"string\"/>"
0415     << endl;
0416   g << "<key id=\"nx\" for=\"node\" attr.name=\"x\" attr.type=\"double\"/>"
0417     << endl;
0418   g << "<key id=\"ny\" for=\"node\" attr.name=\"y\" attr.type=\"double\"/>"
0419     << endl;
0420   g << "<key id=\"nz\" for=\"node\" attr.name=\"z\" attr.type=\"double\"/>"
0421     << endl;
0422   g << "<key id=\"nt\" for=\"node\" attr.name=\"t\" attr.type=\"double\"/>"
0423     << endl;
0424
0425   g << "<key id=\"elabel\" for=\"edge\" attr.name=\"label\" "
0426        "attr.type=\"string\"/>"
0427     << endl;
0428   g << "<key id=\"epl\" for=\"edge\" attr.name=\"plabel\" attr.type=\"int\"/>"
0429     << endl;
0430   g << "<key id=\"epid\" for=\"edge\" attr.name=\"pid\" attr.type=\"int\"/>"
0431     << endl;
0432   g << "<key id=\"estat\" for=\"edge\" attr.name=\"pstat\" attr.type=\"int\"/>"
0433     << endl;
0434
0435   g << "<key id=\"ept\" for=\"edge\" attr.name=\"pT\" attr.type=\"double\"/>"
0436     << endl;
0437   g << "<key id=\"eeta\" for=\"edge\" attr.name=\"eta\" attr.type=\"double\"/>"
0438     << endl;
0439   g << "<key id=\"ephi\" for=\"edge\" attr.name=\"phi\" attr.type=\"double\"/>"
0440     << endl;
0441   g << "<key id=\"ee\" for=\"edge\" attr.name=\"E\" attr.type=\"double\"/>"
0442     << endl;
0443
0444   g << "<graph id=\"G\" edgedefault=\"directed\">" << endl;
0445
0446   node_iterator nIt, nEnd;
0447   int n = 0;
0448
0449   for (nIt = nodes_begin(), nEnd = nodes_end(); nIt != nEnd; ++nIt) {
0450     stringstream stream;
0451
0452     stream << fixed << setprecision(2) << (vMap[*nIt]->x_in().t());
0453
0454     g << "<node id=\"" << n << "\">" << endl;
0455     //g<<"<data key=\"nlabel\">"<<to_string(n)+"("+to_string(vMap[*nIt]->x_in().t())+")"<<"</data>"<<endl;
0456     g << "<data key=\"nlabel\">" << to_string(n) + "(" + stream.str() + ")"
0457       << "</data>" << endl;
0458     g << "<data key=\"nx\">" << vMap[*nIt]->x_in().x() << "</data>" << endl;
0459     g << "<data key=\"ny\">" << vMap[*nIt]->x_in().y() << "</data>" << endl;
0460     g << "<data key=\"nz\">" << vMap[*nIt]->x_in().z() << "</data>" << endl;
0461     g << "<data key=\"nt\">" << vMap[*nIt]->x_in().t() << "</data>" << endl;
0462     g << "</node>" << endl;
0463     n++;
0464   }
0465
0466   edge_iterator eIt, eEnd;
0467   n = 0;
0468   for (eIt = edges_begin(), eEnd = edges_end(); eIt != eEnd; ++eIt) {
0469     g << "<edge id=\"" << n << "\" source=\"" << to_string(eIt->source().id())
0470       << "\" target=\"" << to_string(eIt->target().id()) << "\">" << endl;
0471     g << "<data key=\"elabel\">" << to_string((pMap[*eIt]->pt())) << "</data>"
0472       << endl;
0473     g << "<data key=\"epl\">" << pMap[*eIt]->plabel() << "</data>" << endl;
0474     g << "<data key=\"epid\">" << pMap[*eIt]->pid() << "</data>" << endl;
0475     g << "<data key=\"estat\">" << pMap[*eIt]->pstat() << "</data>" << endl;
0476     g << "<data key=\"ept\">" << pMap[*eIt]->pt() << "</data>" << endl;
0477     g << "<data key=\"eeta\">" << pMap[*eIt]->eta() << "</data>" << endl;
0478     g << "<data key=\"ephi\">" << pMap[*eIt]->phi() << "</data>" << endl;
0479     g << "<data key=\"ee\">" << pMap[*eIt]->e() << "</data>" << endl;
0480     g << "</edge>" << endl;
0481     n++;
0482   }
0483
0484   g << "</graph>" << endl;
0485   g << "</graphml>" << endl;
0486
0487   g.close();
0488
0489 } // end namespace Jetscape
0490 } // namespace Jetscape