HepMC3 event record library
ValidationControl.cc
1 // -*- C++ -*-
2 //
3 // This file is part of HepMC
4 // Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details)
5 //
6 #include "ValidationControl.h"
7 
8 #ifdef SIMPLEEVENTTOOL
9 #include "SimpleEventTool.h"
10 #endif
11 
12 #ifdef PHOTOSPP
13 #include "PhotosValidationTool.h"
14 #endif
15 
16 #ifdef TAUOLAPP
17 #include "TauolaValidationTool.h"
18 #endif
19 
20 #ifdef MCTESTER
21 #include "McTesterValidationTool.h"
22 #endif
23 
24 #ifdef PYTHIA8
25 #include "PythiaValidationTool.h"
26 #endif
27 
28 #include <fstream>
29 #include <cstdio>
30 
31 ValidationControl::ValidationControl():
32 m_events(0),
33 m_momentum_check_events(0),
34 m_momentum_check_threshold(10e-6),
35 m_print_events(0),
36 m_event_counter(0),
37 m_status(-1),
38 m_timer("processing time"),
39 m_has_input_source(0) {
40 }
41 
42 ValidationControl::~ValidationControl() {
43  for (std::vector<ValidationTool *>::iterator t=m_toolchain.begin();t!=m_toolchain.end();++t)
44  delete *t;
45 }
46 
47 void ValidationControl::read_file(const std::string &filename) {
48 
49  // Open config file
50  std::ifstream in(filename.c_str());
51 
52  if(!in.is_open()) {
53  printf("ValidationControl: error reading config file: %s\n",filename.c_str());
54  m_status = -1;
55  return;
56  }
57  else printf("ValidationControl: parsing config file: %s\n",filename.c_str());
58 
59  // Parse config file
60  char buf[256];
61  int line = 0;
62 
63  while(!in.eof()) {
64  PARSING_STATUS status = PARSING_OK;
65  ++line;
66 
67  in >> buf;
68 
69  if( strlen(buf) < 3 || buf[0] == ' ' || buf[0] == '#' ) {
70  in.getline(buf,255);
71  continue;
72  }
73 
74  // Parse event number
75  if( strncmp(buf,"EVENTS",6)==0 ) {
76  in>>m_events;
77  }
78  // Parse input source
79  else if( strncmp(buf,"INPUT",5)==0 ) {
80  in >> buf;
81 
82  if( m_has_input_source ) status = ADDITIONAL_INPUT;
83  else {
84  ValidationTool *input = NULL;
85  // Use tool as input source - currently only one supported tool
86  if( strncmp(buf,"tool",4)==0 ) {
87 #ifdef SIMPLEEVENTTOOL
88  input = new SimpleEventTool();
89 #else
90  status = UNAVAILABLE_TOOL;
91 #endif
92  }
93  else if( strncmp(buf,"pythia8",7)==0) {
94 #ifdef PYTHIA8
95  in >> buf;
96  input = new PythiaValidationTool(buf);
97 #else
98  status = UNAVAILABLE_TOOL;
99 #endif
100  }
101  else status = UNRECOGNIZED_INPUT;
102 
103  if(!status) {
104  m_has_input_source = true;
105  m_toolchain.insert(m_toolchain.begin(),input);
106  }
107  }
108  }
109  // Parse tools used
110  else if( strncmp(buf,"TOOL",3)==0 ) {
111  in >> buf;
112 
113  if ( strncmp(buf,"tauola",6)==0 ) {
114 #ifdef TAUOLAPP
115  m_toolchain.push_back( new TauolaValidationTool() );
116 #else
117  status = UNAVAILABLE_TOOL;
118 #endif
119  }
120  else if( strncmp(buf,"photos",6)==0 ) {
121 #ifdef PHOTOSPP
122  m_toolchain.push_back( new PhotosValidationTool() );
123 #else
124  status = UNAVAILABLE_TOOL;
125 #endif
126  }
127  else if( strncmp(buf,"mctester",8)==0 ) {
128 #ifdef MCTESTER
129  m_toolchain.push_back( new McTesterValidationTool() );
130 #else
131  status = UNAVAILABLE_TOOL;
132 #endif
133  }
134  else status = UNRECOGNIZED_TOOL;
135  }
136  // Parse option
137  else if( strncmp(buf,"SET",3)==0 ) {
138  in >> buf;
139 
140  if ( strncmp(buf,"print_events",12)==0 ) {
141  in >> buf;
142 
143  int events = 0;
144  if( strncmp(buf,"ALL",3)==0 ) events = -1;
145  else events = atoi(buf);
146 
147  print_events(events);
148  }
149  else if( strncmp(buf,"check_momentum",14)==0 ) {
150  in >> buf;
151 
152  int events = 0;
153  if( strncmp(buf,"ALL",3)==0 ) events = -1;
154  else events = atoi(buf);
155 
156  check_momentum_for_events(events);
157  }
158  else status = UNRECOGNIZED_OPTION;
159  }
160  else status = UNRECOGNIZED_COMMAND;
161 
162  // Error checking
163  if(status != PARSING_OK) printf("ValidationControl: config file line %i: ",line);
164 
165  switch(status) {
166  case UNRECOGNIZED_COMMAND: printf("skipping unrecognised command: '%s'\n",buf); break;
167  case UNRECOGNIZED_OPTION: printf("skipping unrecognised option: '%s'\n",buf); break;
168  case UNRECOGNIZED_INPUT: printf("skipping unrecognised input source: '%s'\n",buf); break;
169  case UNRECOGNIZED_TOOL: printf("skipping unrecognised tool: '%s'\n",buf); break;
170  case UNAVAILABLE_TOOL: printf("skipping unavailable tool: '%s'\n",buf); break;
171  case ADDITIONAL_INPUT: printf("skipping additional input source: '%s'\n",buf); break;
172  case CANNOT_OPEN_FILE: printf("skipping input file: '%s'\n",buf); break;
173  default: break;
174  }
175 
176  // Ignore rest of the line
177  in.getline(buf,255);
178  }
179 
180  // Having input source is enough to start validation
181  if(m_has_input_source) m_status = 0;
182  else printf("ValidationControl: no valid input source\n");
183 }
184 
185 bool ValidationControl::new_event() {
186  if( m_status ) return false;
187  if( m_events && ( m_event_counter >= m_events ) ) return false;
188 
189  if(m_event_counter) {
190  if( m_momentum_check_events>0 ) --m_momentum_check_events;
191  if( m_print_events>0 ) --m_print_events;
192  }
193  else m_timer.start();
194 
195  ++m_event_counter;
196 
197  if( m_events ) {
198  if( m_event_counter == 1 ) {
199  printf("ValidationControl: event 1 of %-7i\n",m_events);
200  m_events_print_step = m_events/10;
201  }
202  else if( m_event_counter%m_events_print_step == 0 ) {
203  int elapsed = m_timer.elapsed_time();
204  m_timer.stop();
205  int total = m_timer.total_time();
206  printf("ValidationControl: event %7i (%6.2f%%, %7ims current, %7ims total)\n",m_event_counter,m_event_counter*100./m_events,elapsed,total);
207  m_timer.start();
208  }
209  }
210  else {
211  if( m_event_counter == 1 ) {
212  printf("ValidationControl: event 1\n");
213  m_events_print_step = 1000;
214  }
215  else if( m_event_counter%m_events_print_step == 0 ) {
216  int elapsed = m_timer.elapsed_time();
217  m_timer.stop();
218  int total = m_timer.total_time();
219  printf("ValidationControl: event %7i (%6ims current, %7ims total)\n",m_event_counter,elapsed,total);
220  m_timer.start();
221  }
222  }
223 
224  return true;
225 }
226 
227 void ValidationControl::initialize() {
228  printf("ValidationControl: initializing\n");
229 
230  for (std::vector<ValidationTool *>::iterator tool=m_toolchain.begin();tool!=m_toolchain.end();++tool) (*tool)->initialize();
231 }
232 
233 void ValidationControl::process(GenEvent &hepmc) {
234 
235  m_status = 0;
236 
237  FourVector input_momentum;
238  for (std::vector<ValidationTool *>::iterator tool=m_toolchain.begin();tool!=m_toolchain.end();++tool) {
239 
240  Timer *timer = (*tool)->timer();
241 
242  if(timer) timer->start();
243  m_status = (*tool)->process(hepmc);
244  if(timer) timer->stop();
245 
246  // status != 0 means an error - stop processing current event
247  if(m_status) return;
248 
249  if((*tool)->tool_modifies_event() && m_print_events) {
250  printf("--------------------------------------------------------------\n");
251  printf(" Print event: %s\n",(*tool)->name().c_str());
252  printf("--------------------------------------------------------------\n");
253 
254  HEPMC2CODE( hepmc.print(); )
255  HEPMC3CODE( Print::listing(hepmc,8); )
256  }
257 
258  if((*tool)->tool_modifies_event() && m_momentum_check_events ) {
259  FourVector sum;
260  double delta = 0.0;
261 
262  HEPMC2CODE(
263  for ( GenEvent::particle_const_iterator p = hepmc.particles_begin();
264  p != hepmc.particles_end(); ++p ) {
265  if( (*p)->status() != 1 ) continue;
266  //(*p)->print();
267  FourVector m = (*p)->momentum();
268  sum.setPx( sum.px() + m.px() );
269  sum.setPy( sum.py() + m.py() );
270  sum.setPz( sum.pz() + m.pz() );
271  sum.setE ( sum.e() + m.e() );
272  }
273 
274  double momentum = input_momentum.px() + input_momentum.py() + input_momentum.pz() + input_momentum.e();
275  if( fabs(momentum) > 10e-12 ) {
276  double px = input_momentum.px() - sum.px();
277  double py = input_momentum.py() - sum.py();
278  double pz = input_momentum.pz() - sum.pz();
279  double e = input_momentum.e() - sum.e();
280  delta = sqrt(px*px + py*py + pz*pz + e*e);
281  }
282  )
283 
284  HEPMC3CODE(
285  //vector<GenParticlePtr> results = applyFilter(Selector::STATUS==1,hepmc.particles());
286  for (auto p: hepmc.particles()) if( p->status() != 1 ) continue; else sum += p->momentum();
287  if(!input_momentum.is_zero()) delta = (input_momentum - sum).length();
288  )
289 
290  printf("Momentum sum: %+15.8e %+15.8e %+15.8e %+15.8e (evt: %7i, %s)",sum.px(),sum.py(),sum.pz(),sum.e(),m_event_counter,(*tool)->name().c_str());
291 
292  if( delta < m_momentum_check_threshold ) printf("\n");
293  else printf(" - WARNING! Difference = %+15.8e\n",delta);
294 
295  input_momentum = sum;
296  }
297  }
298 }
299 
300 void ValidationControl::finalize() {
301  printf("ValidationControl: finalizing\n");
302 
303  // Finalize
304  for (std::vector<ValidationTool *>::iterator t=m_toolchain.begin();t!=m_toolchain.end();++t)
305  (*t)->finalize();
306 
307  printf("ValidationControl: printing timers\n");
308 
309  // Print timers
310  for (std::vector<ValidationTool *>::iterator t=m_toolchain.begin();t!=m_toolchain.end();++t)
311  if((*t)->timer()) (*t)->timer()->print();
312 
313 
314  printf("ValidationControl: finished processing:\n");
315 
316  // List tools
317  for (std::vector<ValidationTool *>::iterator t=m_toolchain.begin();t!=m_toolchain.end();++t)
318  printf(" tool: %s\n",(*t)->long_name().c_str());
319 
320 }
void start()
Definition: Timer.h:35
bool is_zero() const
Check if the length of this vertex is zero.
Definition: FourVector.h:154
static void listing(std::ostream &os, const GenEvent &event, unsigned short precision=2)
Print event in listing (HepMC2) format.
Definition: Print.cc:51
void setPy(double pyy)
Set y-component of momentum.
Definition: FourVector.h:89
void stop()
Definition: Timer.h:51
double e() const
Energy component of momentum.
Definition: FourVector.h:97
Stores event-related information.
Definition: GenEvent.h:42
Generic 4-vector.
Definition: FourVector.h:34
double px() const
x-component of momentum
Definition: FourVector.h:82
void setPz(double pzz)
Set z-component of momentum.
Definition: FourVector.h:94
void setPx(double pxx)
Set x-component of momentum.
Definition: FourVector.h:84
double length() const
Magnitude of spatial (x, y, z) 3-vector.
Definition: FourVector.h:110
double pz() const
z-component of momentum
Definition: FourVector.h:92
double py() const
y-component of momentum
Definition: FourVector.h:87
Definition: Timer.h:29
const std::vector< ConstGenParticlePtr > & particles() const
Get list of particles (const)
Definition: GenEvent.cc:40
void setE(double ee)
Set energy component of momentum.
Definition: FourVector.h:99