HepMC3 event record library
HepMC3ViewerFrame.cc
2 #include "HepMC3/ReaderFactory.h"
3 #include "HepMC3/GenEvent.h"
4 #include "HepMC3/GenParticle.h"
5 #include "HepMC3/GenVertex.h"
6 #include "HepMC3ViewerFrame.h"
7 
8 /* Older graphviz versions can have conflictiong declarations of memcmp/strcmp function
9  * This can break compilation with -pedantic. Uncomenting line below can fix it.
10  */
11 // #define _PACKAGE_ast 1
12 
13 #include <graphviz/gvc.h>
14 #define CONSERVATION_TOLERANCE 1e-5
15 
16 static char* create_image_from_dot(char* m_buffer)
17 {
18  GVC_t * gvc=gvContext();
19  Agraph_t *g= agmemread(m_buffer);
20  gvLayout(gvc,g,"dot");
21 
22  int err;
23  char *data;
24  unsigned int length;
25 
26  if (!g)
27  return NULL;
28  err = gvRenderData(gvc, g, "png", &data, &length);
29  if (err)
30  return NULL;
31  data = (char*)realloc(data, length + 1);
32  delete g;
33  delete gvc;
34  return data;
35 }
36 
37 static bool show_as_parton(HepMC3::ConstGenParticlePtr p )
38 {
39  const int pd=std::abs(p->pid());
40  bool parton=false;
41 
42  if (pd==81||pd==82||pd<25) parton=true;
43  if (
44  (pd/1000==1||pd/1000==2||pd/1000==3||pd/1000==4||pd/1000==5)
45  &&(pd%1000/100==1||pd%1000/100==2||pd%1000/100==3||pd%1000/100==4)
46  &&(pd%100==1||pd%100==3)
47  )
48  parton=true;
49  if (p->status()==4) parton=true;
50  return parton;
51 }
52 
53 static char* write_event_to_dot(char* used_cursor,const HepMC3::GenEvent &evt,int used_style=1)
54 {
55  used_cursor += sprintf(used_cursor, "digraph graphname%d {\n",evt.event_number());
56  used_cursor += sprintf(used_cursor, "v0[label=\"Machine\"];\n");
57  for(auto v: evt.vertices() )
58  {
59  if (used_style!=0)
60  {
61  if (used_style==1) //paint decay and fragmentation vertices in green
62  {
63  if (v->status()==2) used_cursor += sprintf(used_cursor, "node [color=\"green\"];\n");
64  else used_cursor += sprintf(used_cursor, "node [color=\"black\"];\n");
65  }
66  }
69  double energy=0;
70  for(auto p1: v->particles_in() ) {
71  in+=p1->momentum();
72  energy+=std::abs(p1->momentum().e());
73  }
74  for(auto p2: v->particles_out() ) {
75  out+=p2->momentum();
76  energy+=std::abs(p2->momentum().e());
77  }
78  HepMC3::FourVector momviolation(0,0,0,0);
79  momviolation+=in;
80  momviolation-=out;
81  double energyviolation=std::sqrt(momviolation.length2() +momviolation.e()*momviolation.e() );
82  bool violation=false;
83  if (energyviolation>CONSERVATION_TOLERANCE*energy) violation=true;
84 
85  if(violation)
86  {
87  used_cursor += sprintf(used_cursor, "node [shape=rectangle];\n");
88  used_cursor += sprintf(used_cursor, "v%d [label=\"%d\nd=%4.2f\"];\n", -v->id(),v->id(),energyviolation);
89  }
90  else
91  {
92  used_cursor += sprintf(used_cursor, "node [shape=ellipse];\n");
93  used_cursor += sprintf(used_cursor, "v%d[label=\"%d\"];\n", -v->id(),v->id());
94  }
95 
96  used_cursor += sprintf(used_cursor, "node [shape=ellipse];\n");
97  }
98  for(auto p: evt.beams() )
99  {
100  if (!p->end_vertex()) continue;
101  used_cursor += sprintf(used_cursor, "node [shape=point];\n");
102  used_cursor += sprintf(used_cursor, "v0 -> v%d [label=\"%d(%d)\"];\n", -p->end_vertex()->id(),p->id(),p->pid());
103  }
104 
105  for(auto v: evt.vertices() )
106  {
107 
108  for(auto p: v->particles_out() )
109  {
110  {
111  if (used_style!=0)
112  {
113  if (used_style==1) //paint suspected partons and 81/82 in red
114  {
115  if (show_as_parton(p)&&p->status()!=1) used_cursor += sprintf(used_cursor, "edge [color=\"red\"];\n");
116  else used_cursor +=sprintf(used_cursor, "edge [color=\"black\"];\n");
117  }
118  }
119  if (!p->end_vertex())
120  {
121  used_cursor += sprintf(used_cursor, "node [shape=point];\n");
122  used_cursor += sprintf(used_cursor, "v%d -> o%d [label=\"%d(%d)\"];\n", -v->id(),p->id(),p->id(),p->pid());
123  continue;
124  }
125  else
126  used_cursor += sprintf(used_cursor, "v%d -> v%d [label=\"%d(%d)\"];\n", -v->id(),-p->end_vertex()->id(),p->id(),p->pid());
127  }
128  }
129  }
130  used_cursor += sprintf(used_cursor, "labelloc=\"t\";\nlabel=\"Event %d; Vertices %lu; Particles %lu;\";\n", evt.event_number(), evt.vertices().size(), evt.particles().size());
131  used_cursor += sprintf(used_cursor,"}\n\n");
132 
133  return used_cursor;
134 }
135 
136 
138 {
139  char* m_buffer = new char[m_char_buffer_size]();
140  char* m_cursor=m_buffer;
141  m_cursor=write_event_to_dot(m_cursor,*(fCurrentEvent));
142  char *buf=create_image_from_dot(m_buffer);
143  fEmbEventImageCanvas->MapSubwindows();
144 
145  if(!fEventImageCanvas) fEventImageCanvas=new TCanvas("fEmbEventImageCanvas","fEmbEventImageCanvas",1024,768);
146 
147  fEventImageCanvas->cd();
148  fEventImageCanvas->Clear();
149  double d=0.60;
150 
151  fGraphImage = TImage::Create();
152  fGraphImage->SetName("Event");
153  fGraphImage->SetImageBuffer(&buf, TImage::kPng);
154 
155  fGraphImage->SetConstRatio(kFALSE);
156 
157  TPad *p1 = new TPad("i1", "i1", 0.05, 0.05, 0.05+d*fGraphImage->GetWidth()/fGraphImage->GetHeight(), 0.95);
158  p1->Draw();
159  p1->cd();
160 
161  fGraphImage->Draw("xxx");
162  delete [] m_buffer;
163  gPad->Update();
164  DoAnalysis();
165 }
166 
168 {
169  fEmbAnalysisCanvas->MapSubwindows();
170  fAnalysisCanvas->cd();
171  fAnalysisCanvas->Clear();
172  for (auto h: fAnalysisH) h.second->Delete();
173  fAnalysisH.clear();
174 
175  /* */
176  TH1S* particles1= new TH1S();
177  fAnalysisH["particles1"]=particles1;
178  particles1->SetTitle("Flavour: all particles; PDG ID; Number of particles");
179  particles1->SetFillColor(kBlue);
180  for(auto p: fCurrentEvent->particles() )
181  particles1->Fill((std::to_string(p->pid())).c_str(),1.0);
182  particles1->LabelsOption(">","X");
183  /* */
184  TH1S* particles2= new TH1S();
185  fAnalysisH["particles2"]=particles2;
186  particles2->SetTitle("Flavour: particles with status 1; PDG ID; Number of particles");
187  particles2->SetFillColor(kBlue);
188  for(auto p: fCurrentEvent->particles() )
189  if(p->status()==1) particles2->Fill((std::to_string(p->pid())).c_str(),1.0);
190  particles2->LabelsOption(">","X");
191  /* */
192  std::vector<double> masses;
193  for(auto p: fCurrentEvent->particles() )
194  if(show_as_parton(p)) masses.push_back(p->momentum().m());
195  TH1D* particles3= new TH1D("particles3","Mass: parton particles; Mass, GeV; Number of particles",masses.size(),
196  0,*std::max_element(masses.begin(),masses.end()));
197  fAnalysisH["particles3"]=particles3;
198  particles3->SetFillColor(kBlue);
199  for(auto m: masses) particles3->Fill(m);
200 
201 
202  fAnalysisCanvas->cd();
203  TPad *p1 = new TPad("i1", "i1", 0.00, 0.75, 1.0, 1.0);
204  p1->Draw();
205  p1->cd();
206  particles1->Draw();
207  fAnalysisCanvas->cd();
208  TPad *p2 = new TPad("i2", "i2", 0.00, 0.50, 1.0, 0.75);
209  p2->Draw();
210  p2->cd();
211  particles2->Draw();
212  fAnalysisCanvas->cd();
213  TPad *p3 = new TPad("i3", "i3", 0.00, 0.25, 1.0, 0.50);
214  p3->Draw();
215  p3->cd();
216  particles3->Draw();
217 
218  gPad->Update();
219 }
220 
222 {
223  fEventsCache.clear();
224  fCurrentEvent=nullptr;
225 }
226 
228 {
229  auto pos=find(fEventsCache.begin(),fEventsCache.end(),fCurrentEvent);
230  if (pos==fEventsCache.begin()) return;
231  pos--;
232  fCurrentEvent=*(pos);
233  if (pos==fEventsCache.end()) printf("This event was not found in the cache. Cache size is %zu\n",fEventsCache.size());
234  DrawEvent();
235 }
236 
237 void HepMC3ViewerFrame::ReadFile(const char* a) {
239 }
240 
242 {
243  if (fCurrentEvent==nullptr||fEventsCache.back()==fCurrentEvent)
244  {
245  HepMC3::GenEvent* evt1=new HepMC3::GenEvent(HepMC3::Units::GEV,HepMC3::Units::MM);
246  bool ok=fReader->read_event(*(evt1));
247  ok=(ok&&!fReader->failed());
248  if (ok)
249  {
250  fEventsCache.push_back(evt1);
251  fCurrentEvent=evt1;
252  }
253  else return;
254  }
255  else
256  {
257  auto pos=find(fEventsCache.begin(),fEventsCache.end(),fCurrentEvent);
258  pos++;
259  fCurrentEvent=*(pos);
260  }
261  DrawEvent();
262 }
264 {
265  static const char *FileType[] = {"All", "*.*","HepMC", "*.hepmc*","LHEF", "*.lhe*","ROOT", "*.root", 0, 0 };
266  static TString dir("./");
267  TGFileInfo fi;
268  fi.fFileTypes = FileType;
269  fi.fIniDir = StrDup(dir);
270  new TGFileDialog(gClient->GetRoot(), this, kFDOpen, &fi);
271  if (fReader) fReader->close();
272  fReader=HepMC3::deduce_reader(fi.fFilename);
273 }
274 
275 HepMC3ViewerFrame::HepMC3ViewerFrame(const TGWindow *p, UInt_t w, UInt_t h) :
276  TGMainFrame(p, w, h)
277 {
278  fMainFrame = new TGCompositeFrame(this, 1350, 500, kHorizontalFrame|kFixedWidth);
279  fButtonFrame = new TGCompositeFrame(fMainFrame, 150, 200, kFixedWidth);
280 
281  fEmbEventImageCanvas =new TRootEmbeddedCanvas("MainCanvaslegent", fMainFrame, 850, 500);
282 
283  fEmbAnalysisCanvas =new TRootEmbeddedCanvas("EmbAnalysisCanvaslegend", fMainFrame, 350, 500);
284 
285 
286  fMainFrame->AddFrame(fEmbEventImageCanvas,new TGLayoutHints(kLHintsTop | kLHintsExpandX| kLHintsExpandY, 1, 1, 2, 2));
287  fMainFrame->AddFrame(fEmbAnalysisCanvas,new TGLayoutHints(kLHintsTop | kFixedWidth| kLHintsExpandY, 1, 1, 2, 2));
288  fMainFrame->AddFrame(fButtonFrame,new TGLayoutHints(kLHintsTop, 1, 1, 2, 2));
289 
290 
291  fChooseInput = new TGTextButton(fButtonFrame, "&Choose input");
292  fChooseInput->Connect("Clicked()", "HepMC3ViewerFrame", this, "ChooseInput()");
293  fChooseInput->SetToolTipText("Click to choose file");
294  fButtonFrame->AddFrame(fChooseInput, new TGLayoutHints(kLHintsTop | kLHintsExpandX, 1, 1, 2, 2));
295 
296 
297 
298  fNextEvent = new TGTextButton(fButtonFrame, "&Next event");
299  fNextEvent->Connect("Clicked()", "HepMC3ViewerFrame", this, "NextEvent()");
300  fNextEvent->SetToolTipText("Click to display next event");
301  fButtonFrame->AddFrame(fNextEvent, new TGLayoutHints(kLHintsExpandX|kLHintsLeft, 1, 1, 2, 2));
302 
303 
304  fPreviousEvent = new TGTextButton(fButtonFrame, "&Previous event");
305  fPreviousEvent->Connect("Clicked()", "HepMC3ViewerFrame", this, "PreviousEvent()");
306  fPreviousEvent->SetToolTipText("Click to display previous event");
307  fButtonFrame->AddFrame(fPreviousEvent, new TGLayoutHints( kLHintsExpandX|kLHintsLeft, 1, 1, 2, 2));
308 
309 
310  fClearEventCache = new TGTextButton(fButtonFrame, "&Clear event cache");
311  fClearEventCache->Connect("Clicked()", "HepMC3ViewerFrame", this, "ClearEventCache()");
312  fClearEventCache->SetToolTipText("Click to clear event cache ");
313  fButtonFrame->AddFrame(fClearEventCache, new TGLayoutHints( kLHintsExpandX|kLHintsLeft, 1, 1, 2, 2));
314 
315  fExit = new TGTextButton(fButtonFrame, "&Exit ","gApplication->Terminate(0)");
316  fExit->SetToolTipText("Click to exit");
317  fButtonFrame->AddFrame(fExit, new TGLayoutHints( kLHintsExpandX|kLHintsLeft,1,1,2,2));
318 
319  AddFrame(fMainFrame, new TGLayoutHints(kLHintsTop |kLHintsExpandX| kLHintsExpandY, 1, 1, 2, 2));
320 
321  SetWindowName("Event viewer");
322  MapSubwindows();
323  Resize(GetDefaultSize());
324  MapWindow();
325 
326  fReader=nullptr;
328  fAnalysisCanvas=fEmbAnalysisCanvas->GetCanvas();
329  fCurrentEvent=nullptr;
330  fGraphImage = TImage::Create();
331 }
333 {
334  fMainFrame->Cleanup();
335  fReader->close();
336  Cleanup();
337 }
int event_number() const
Get event number.
Definition: GenEvent.h:136
TGCompositeFrame * fButtonFrame
Button frame.
std::map< std::string, TH1 * > fAnalysisH
Analysis histograms.
TGTextButton * fClearEventCache
Button.
Definition of class GenParticle.
const std::vector< ConstGenVertexPtr > & vertices() const
Get list of vertices (const)
Definition: GenEvent.cc:43
TCanvas * fEventImageCanvas
Event canvas.
void DrawEvent()
Draw evemt.
Definition of class GenVertex.
Definition of class ReaderRootTree.
std::shared_ptr< Reader > deduce_reader(std::istream &stream)
This function will deduce the type of input stream based on its content and will return appropriate R...
TGTextButton * fNextEvent
Button.
TCanvas * fAnalysisCanvas
Analysis canvas.
HepMC3::GenEvent * fCurrentEvent
Event.
TRootEmbeddedCanvas * fEmbEventImageCanvas
Event canvas.
std::shared_ptr< HepMC3::Reader > fReader
Reader.
Stores event-related information.
Definition: GenEvent.h:41
Generic 4-vector.
Definition: FourVector.h:36
std::vector< HepMC3::GenEvent * > fEventsCache
Cache of events.
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
Definition: GenEvent.cc:39
void ClearEventCache()
slot
TGTextButton * fExit
Button.
virtual ~HepMC3ViewerFrame()
Destructor.
TGTextButton * fChooseInput
Button.
static const size_t m_char_buffer_size
Size of writer buffer.
TImage * fGraphImage
Image passed from graphviz.
TGCompositeFrame * fMainFrame
Main frame.
TRootEmbeddedCanvas * fEmbAnalysisCanvas
Analysis canvas.
Definition of class GenEvent.
HepMC3ViewerFrame(const TGWindow *p, UInt_t w, UInt_t h)
Constructor.
void DoAnalysis()
Do analysis.
void ReadFile(const char *a)
Open file.
TGTextButton * fPreviousEvent
Button.
std::vector< ConstGenParticlePtr > beams() const
Vector of beam particles.
Definition: GenEvent.cc:416
Feature< Feature_type > abs(const Feature< Feature_type > &input)
Obtain the absolute value of a Feature. This works as you&#39;d expect. If foo is a valid Feature...
Definition: Feature.h:323