PYTHIA  8.311
Visualisation.h
1 // Visualisation.h is a part of the PYTHIA event generator.
2 // Copyright (C) 2024 Torbjorn Sjostrand.
3 // PYTHIA is licenced under the GNU GPL v2 or later, see COPYING for details.
4 // Please respect the MCnet Guidelines, see GUIDELINES for details.
5 
6 // The following functions analyze a scattering event and save the event in
7 // an output format that can be converted into a postscript figure using the
8 // "graphviz" program.
9 // Author: Nadine Fischer (primary), Stefan Prestel (porting)
10 
11 #ifndef Pythia8_Visualisation_H
12 #define Pythia8_Visualisation_H
13 
14 namespace Pythia8 {
15 
16 //==========================================================================
17 
18 // The functions below assist the main function, printEvent, in producing
19 // an event visualisation, by converting status codes and defining arrows.
20 
21 string py_status(int stAbs) {
22  string status = "";
23  if (stAbs > 20 && stAbs < 30) status = "hardProcess";
24  else if (stAbs > 30 && stAbs < 40) status = "MPI";
25  else if (stAbs > 40 && stAbs < 50) status = "ISR";
26  else if (stAbs > 50 && stAbs < 60) status = "FSR";
27  else if (stAbs > 60 && stAbs < 70) status = "beamRemnants";
28  else if (stAbs > 70 && stAbs < 80) status = "hadronizationPrep";
29  else if (stAbs > 80 && stAbs < 90) status = "hadronization";
30  else if (stAbs > 90 && stAbs < 110) status = "decays";
31  else status = "default";
32  return status;
33 }
34 
35 //--------------------------------------------------------------------------
36 
37 void makeArrow(map< pair<string,string>, string >* arrows,
38  string identParent, string identChild) {
39  pair<string,string> key = make_pair(identParent,identChild);
40  string value = " " + identParent + " -> " + identChild
41  + " [weight=2,label=\" \"];";
42  arrows->insert( pair< pair<string,string>, string>(key, value) );
43 }
44 
45 //--------------------------------------------------------------------------
46 
47 // The main visualisation function. Note that this is really only a schematic
48 // visualisation. In particular, the space-time structure of the partonic
49 // evolution and the space-time position of string breakups are only handled
50 // schematically.
51 
52 void printEvent(Event& evt, string fileName = "event") {
53 
54  bool simplifyHadronization = true;
55  bool addLegend = true;
56  map<string, pair<string,string> > colMap;
57  colMap["default"] = make_pair("white","black");
58  colMap["hardProcess"] = make_pair("red","black");
59  colMap["MPI"] = make_pair("lightsalmon","black");
60  colMap["ISR"] = make_pair("lightseagreen","black");
61  colMap["FSR"] = make_pair("limegreen","black");
62  colMap["beamRemnants"] = make_pair("mediumpurple","black");
63  colMap["hadronizationPrep"] = make_pair("blue","black");
64  colMap["hadronization"] = make_pair("blue","black");
65  colMap["decays"] = make_pair("lightskyblue","black");
66 
67  map<string,string> blobs;
68  map< pair<string,string>, string > arrows;
69  vector< vector<int> > hadronGroups;
70  vector< vector<int> > hadronParents;
71 
72  for (int i=1; i<(int)evt.size(); i++) {
73  // Identifier of that particle.
74  string ident = "F" + std::to_string(10000+i);
75  // Name that will appear in graph.
76  string label = std::to_string(evt[i].id()) + " (" + evt[i].name() + ")";
77  // Find particle group for colors.
78  string status = py_status(evt[i].statusAbs());
79  // Skip hadrons and decay products for simplified output.
80  if (simplifyHadronization &&
81  (status == "decays" || status == "hadronization") ) continue;
82  // Special treatment of hadronization particles for simplified output.
83  bool checkDaughters = simplifyHadronization;
84  if (status != "hadronizationPrep" && status != "beamRemnants")
85  checkDaughters = false;
86  // Check that daughters are are part of hadronization
87  if (checkDaughters) {
88  vector<int> daus = evt[i].daughterList();
89  for (int j=0; j<(int)daus.size(); j++)
90  if (py_status(evt[daus[j]].statusAbs()) != "hadronization")
91  checkDaughters = false;
92  }
93  if (checkDaughters) {
94  vector<int> daus = evt[i].daughterList();
95  // Check if other particles in preparation has same daughter list.
96  bool foundSameDaus = false;
97  for (int j=0; j<(int)hadronGroups.size(); j++) {
98  if (daus.size() == hadronGroups[j].size()) {
99  foundSameDaus = true;
100  for (int k=0; k<(int)hadronGroups[j].size(); k++)
101  if (daus[k] != hadronGroups[j][k]) foundSameDaus = false;
102  if (foundSameDaus) {
103  hadronParents[j].push_back(i);
104  break;
105  }
106  }
107  }
108  if (!foundSameDaus) {
109  hadronGroups.push_back(daus);
110  vector<int> parents; parents.push_back(i);
111  hadronParents.push_back(parents);
112  }
113  if (status == "hadronizationPrep") continue;
114  }
115  // Setup the graph for the particle.
116  pair<string,string> colors = colMap[status];
117  string fillcolor = colors.first, fontcolor = colors.second;
118  blobs[ident] = " " + ident + " [shape=box,style=filled,fillcolor=\""
119  + fillcolor + "\",fontcolor=\"" + fontcolor + "\",label=\""
120  + label + "\"];";
121  // Setup arrow to mother(s).
122  int mot1 = evt[i].mother1(), mot2 = evt[i].mother2();
123  if ( i > 3 && (mot1 == 0 || mot2 == 0) )
124  makeArrow(&arrows, "F"+std::to_string(10000+max(mot1,mot2)), ident);
125  // Setup arrow to daughter(s).
126  if (!checkDaughters) {
127  vector<int> daus = evt[i].daughterList();
128  for (int j=0; j<(int)daus.size(); j++)
129  makeArrow(&arrows, ident, "F"+std::to_string(10000+daus[j]));
130  }
131  }
132 
133  // Add the hadron groups for simplified output.
134  map< pair<string,string>, string > arrowsSav = arrows;
135  for (int i=0; i<(int)hadronGroups.size(); i++) {
136  // Identifier of that group.
137  string ident = "G" + std::to_string(10000+i);
138  pair<string,string> colors = colMap["hadronization"];
139  string fillcolor = colors.first, fontcolor = colors.second;
140  string line = " " + ident + " [shape=none,\n label = <<"
141  "table border=\"0\" cellspacing=\"0\">\n";
142  for (int j=0; j<(int)hadronGroups[i].size(); j++) {
143  // Name that will appear in graph.
144  string label = std::to_string(evt[hadronGroups[i][j]].id()) + " ("
145  + evt[hadronGroups[i][j]].name() + ")";
146  line += ( " <tr><td port=\"port" + std::to_string(j)
147  + "\" border=\"1\" bgcolor=\"" + fillcolor + "\"><font color=\""
148  + fontcolor + "\">" + label + "</font></td></tr>\n" );
149  }
150  line += " </table>> ];";
151  // Add the group to the graph.
152  blobs[ident] = line;
153  // Add an arrow from each parent to the group.
154  for (int j=0; j<(int)hadronParents[i].size(); j++) {
155  // Identifier of that parent.
156  string identParent = "F"+std::to_string(10000+hadronParents[i][j]);
157  // List of particles to be erased.
158  vector<string> toErase;
159  toErase.push_back(identParent);
160  // Check if parent is beam remnant.
161  bool parentIsBR = (py_status(evt[hadronParents[i][j]].statusAbs()) ==
162  "beamRemnants");
163  if (parentIsBR) {
164  makeArrow(&arrows, identParent, ident);
165  } else {
166  int nrGP1 = evt[hadronParents[i][j]].mother1();
167  int nrGP2 = evt[hadronParents[i][j]].mother2();
168  if (nrGP1 > 0) {
169  // Trace back one more generation if double hadronization prep.
170  if (py_status(evt[nrGP1].statusAbs()) == "hadronizationPrep") {
171  toErase.push_back("F"+std::to_string(10000+nrGP1));
172  int nrGGP1 = evt[nrGP1].mother1();
173  int nrGGP2 = evt[nrGP1].mother2();
174  if (nrGGP1 > 0)
175  makeArrow(&arrows, "F"+std::to_string(10000+nrGGP1), ident);
176  if (nrGGP2 > 0 && nrGGP2 != nrGGP1)
177  makeArrow(&arrows, "F"+std::to_string(10000+nrGGP2), ident);
178  } else makeArrow(&arrows, "F"+std::to_string(10000+nrGP1), ident);
179  }
180  if (nrGP2 > 0 && nrGP2 != nrGP1) {
181  // Trace back one more generation if double hadronization prep.
182  if (py_status(evt[nrGP2].statusAbs()) == "hadronizationPrep") {
183  toErase.push_back("F"+std::to_string(10000+nrGP2));
184  int nrGGP1 = evt[nrGP2].mother1();
185  int nrGGP2 = evt[nrGP2].mother2();
186  if (nrGGP1 > 0)
187  makeArrow(&arrows, "F"+std::to_string(10000+nrGGP1), ident);
188  if (nrGGP2 > 0 && nrGGP2 != nrGGP1)
189  makeArrow(&arrows, "F"+std::to_string(10000+nrGGP2), ident);
190  } else makeArrow(&arrows, "F"+std::to_string(10000+nrGP2), ident);
191  }
192  // Erase any parents that might be left in the graph.
193  for (int iToE=0; iToE<(int)toErase.size(); iToE++)
194  if (blobs.find(toErase[iToE]) != blobs.end())
195  blobs.erase(toErase[iToE]);
196  for (map< pair<string,string>, string >::iterator k=arrowsSav.begin();
197  k!=arrowsSav.end(); k++) {
198  for (int iToE=0; iToE<(int)toErase.size(); iToE++) {
199  if (k->first.second == toErase[iToE])
200  arrows.erase(k->first);
201  }
202  }
203  }
204  }
205  }
206 
207  // Write output.
208  ofstream outfile;
209  outfile.open((char*)(fileName+".dot").c_str());
210  outfile << "digraph \"event\" {" << endl
211  << " rankdir=LR;" << endl;
212  for (map<string,string>::iterator iBlob=blobs.begin(); iBlob!=blobs.end();
213  iBlob++) outfile << iBlob->second << endl;
214  for (map< pair<string,string>, string >::iterator iArrow=arrows.begin();
215  iArrow!=arrows.end(); iArrow++) outfile << iArrow->second << endl;
216  // Add a legend, skip default.
217  if (addLegend) {
218  outfile << " { rank = source;" << endl
219  << " Legend [shape=none, margin=0, label=<<table border=\"0\""
220  << " cellspacing=\"0\">" << endl
221  << " <tr><td port=\"0\" border=\"1\"><b>Legend</b></td></tr>"
222  << endl;
223  int count = 1;
224  for (map<string, pair<string,string> >::iterator iLeg=colMap.begin();
225  iLeg!=colMap.end(); iLeg++) {
226  if (iLeg->first == "default") continue;
227  if (iLeg->first == "hadronizationPrep") continue;
228  if (simplifyHadronization && iLeg->first == "decays") continue;
229  string fillcolor = iLeg->second.first;
230  string fontcolor = iLeg->second.second;
231  outfile << " <tr><td port=\"port" << count << "\" border=\"1\" "
232  << "bgcolor=\"" << fillcolor << "\"><font color=\"" << fontcolor
233  << "\">" << iLeg->first << "</font></td></tr>" << endl;
234  count++;
235  }
236  outfile << " </table>" << endl << " >];" << endl << " }" << endl;
237  }
238  outfile << "}" << endl;
239  outfile.close();
240 
241  cout << "\n\nPrinted one event to output file " << fileName + ".dot\n";
242  if (system(nullptr)) {
243  if (system("which dot > /dev/null 2>&1") == 0) {
244  cout << "Producing .ps figure by using the 'dot' command." << endl;
245  string command = "dot -Tps " + fileName + ".dot -o " + fileName+".ps";
246  if (system(command.c_str()) == 0)
247  cout << "Stored event visualization in file " << fileName+".ps"
248  << endl;
249  else
250  cout << "Failed to store event visualization in file." << endl;
251  }
252  } else {
253  cout << "You can now produce a .ps figure by using the 'dot' command:\n\n"
254  << "dot -Tps " << fileName << ".dot -o " << fileName << ".ps" << "\n\n";
255  cout << "Note: 'dot' is part of the 'graphviz' package.\n"
256  << "You might want to install this package to produce the .ps event"
257  << endl << endl;
258  }
259 }
260 
261 //==========================================================================
262 
263 } // end namespace Pythia8
264 
265 #endif // end Pythia8_Visualisation_H
The Event class holds all info on the generated event.
Definition: Event.h:453
void printEvent(Event &evt, string fileName="event")
Definition: Visualisation.h:52
string py_status(int stAbs)
Definition: Visualisation.h:21
int size() const
Event record size.
Definition: Event.h:504
Header for classes to set beam momentum and interaction vertex spread.
Definition: Analysis.h:20
vector< int > daughterList(int i) const
Definition: Event.h:593