Fork me on GitHub

source: git/external/HepMC3/ReaderAsciiHepMC2.cc@ 95a917c

Last change on this file since 95a917c was 95a917c, checked in by Pavel Demin <pavel.demin@…>, 3 years ago

add HepMC3 library
  • Property mode set to 100644
File size: 25.5 KB
RevLine 
[95a917c]1// -*- C++ -*-
2//
3// This file is part of HepMC
4// Copyright (C) 2014-2019 The HepMC collaboration (see AUTHORS for details)
5//
6/**
7 * @file ReaderAsciiHepMC2.cc
8 * @brief Implementation of \b class ReaderAsciiHepMC2
9 *
10 */
11#include "HepMC3/ReaderAsciiHepMC2.h"
12
13#include "HepMC3/GenEvent.h"
14#include "HepMC3/GenVertex.h"
15#include "HepMC3/GenParticle.h"
16#include "HepMC3/GenHeavyIon.h"
17#include "HepMC3/GenPdfInfo.h"
18#include "HepMC3/Setup.h"
19#include <cstring>
20#include <cstdlib>
21namespace HepMC3 {
22
23ReaderAsciiHepMC2::ReaderAsciiHepMC2(const std::string& filename):
24 m_file(filename), m_stream(0), m_isstream(false) {
25 if( !m_file.is_open() ) {
26 HEPMC3_ERROR( "ReaderAsciiHepMC2: could not open input file: "<<filename )
27 }
28 set_run_info(std::make_shared<GenRunInfo>());
29 m_event_ghost= new GenEvent();
30}
31// Ctor for reading from stdin
32ReaderAsciiHepMC2::ReaderAsciiHepMC2(std::istream & stream)
33 : m_stream(&stream), m_isstream(true)
34{
35 if( !m_stream->good() ) {
36 HEPMC3_ERROR( "ReaderAsciiHepMC2: could not open input stream " )
37 }
38 set_run_info(std::make_shared<GenRunInfo>());
39 m_event_ghost= new GenEvent();
40}
41
42ReaderAsciiHepMC2::~ReaderAsciiHepMC2() { if (m_event_ghost) { m_event_ghost->clear(); delete m_event_ghost; m_event_ghost=nullptr; } if (!m_isstream) close(); }
43
44bool ReaderAsciiHepMC2::skip(const int n)
45{
46 const size_t max_buffer_size=512*512;
47 char buf[max_buffer_size];
48 int nn=n;
49 while(!failed()) {
50 char peek;
51 if ( (!m_file.is_open()) && (!m_isstream) ) return false;
52 m_isstream ? peek = m_stream->peek() : peek = m_file.peek();
53 if( peek=='E' ) nn--;
54 if (nn<0) return true;
55 m_isstream ? m_stream->getline(buf,max_buffer_size) : m_file.getline(buf,max_buffer_size);
56 }
57 return true;
58}
59
60bool ReaderAsciiHepMC2::read_event(GenEvent &evt) {
61 if ( (!m_file.is_open()) && (!m_isstream) ) return false;
62
63 char peek;
64 const size_t max_buffer_size=512*512;
65 const size_t max_weights_size=256;
66 char buf[max_buffer_size];
67 bool parsed_event_header = false;
68 bool is_parsing_successful = true;
69 int parsing_result = 0;
70 unsigned int vertices_count = 0;
71 unsigned int current_vertex_particles_count = 0;
72 unsigned int current_vertex_particles_parsed= 0;
73
74 evt.clear();
75 evt.set_run_info(run_info());
76
77 // Empty cache
78 m_vertex_cache.clear();
79 m_vertex_barcodes.clear();
80
81 m_particle_cache.clear();
82 m_end_vertex_barcodes.clear();
83 m_particle_cache_ghost.clear();
84 //
85 // Parse event, vertex and particle information
86 //
87 while(!failed()) {
88 m_isstream ? m_stream->getline(buf,max_buffer_size) : m_file.getline(buf,max_buffer_size);
89 if( strlen(buf) == 0 ) continue;
90 // Check for IO_GenEvent header/footer
91 if( strncmp(buf,"HepMC",5) == 0 ) {
92 if( strncmp(buf,"HepMC::Version",14) != 0 && strncmp(buf,"HepMC::IO_GenEvent",18)!=0 )
93 {
94 HEPMC3_WARNING( "ReaderAsciiHepMC2: found unsupported expression in header. Will close the input." )
95 std::cout<<buf<<std::endl;
96 m_isstream ? m_stream->clear(std::ios::eofbit) : m_file.clear(std::ios::eofbit);
97 }
98 if(parsed_event_header) {
99 is_parsing_successful = true;
100 break;
101 }
102 continue;
103 }
104 switch(buf[0]) {
105 case 'E':
106 parsing_result = parse_event_information(evt,buf);
107 if(parsing_result<0) {
108 is_parsing_successful = false;
109 HEPMC3_ERROR( "ReaderAsciiHepMC2: HEPMC3_ERROR parsing event information" )
110 }
111 else {
112 vertices_count = parsing_result;
113 m_vertex_cache.reserve(vertices_count);
114 m_particle_cache.reserve(vertices_count*3);
115 m_vertex_barcodes.reserve(vertices_count);
116 m_end_vertex_barcodes.reserve(vertices_count*3);
117 is_parsing_successful = true;
118 }
119 parsed_event_header = true;
120 break;
121 case 'V':
122 // If starting new vertex: verify if previous was fully parsed
123
124 /** @bug HepMC2 files produced with Pythia8 are known to have wrong
125 information about number of particles in vertex. Hence '<' sign */
126 if(current_vertex_particles_parsed < current_vertex_particles_count) {
127 is_parsing_successful = false;
128 break;
129 }
130 current_vertex_particles_parsed = 0;
131
132 parsing_result = parse_vertex_information(buf);
133
134 if(parsing_result<0) {
135 is_parsing_successful = false;
136 HEPMC3_ERROR( "ReaderAsciiHepMC2: HEPMC3_ERROR parsing vertex information" )
137 }
138 else {
139 current_vertex_particles_count = parsing_result;
140 is_parsing_successful = true;
141 }
142 break;
143 case 'P':
144
145 parsing_result = parse_particle_information(buf);
146
147 if(parsing_result<0) {
148 is_parsing_successful = false;
149 HEPMC3_ERROR( "ReaderAsciiHepMC2: HEPMC3_ERROR parsing particle information" )
150 }
151 else {
152 ++current_vertex_particles_parsed;
153 is_parsing_successful = true;
154 }
155 break;
156 case 'U':
157 is_parsing_successful = parse_units(evt,buf);
158 break;
159 case 'F':
160 is_parsing_successful = parse_pdf_info(evt,buf);
161 break;
162 case 'H':
163 is_parsing_successful = parse_heavy_ion(evt,buf);
164 break;
165 case 'N':
166 is_parsing_successful = parse_weight_names(buf);
167 break;
168 case 'C':
169 is_parsing_successful = parse_xs_info(evt,buf);
170 break;
171 default:
172 HEPMC3_WARNING( "ReaderAsciiHepMC2: skipping unrecognised prefix: " << buf[0] )
173 is_parsing_successful = true;
174 break;
175 }
176
177 if( !is_parsing_successful ) break;
178
179 // Check for next event
180 m_isstream ? peek = m_stream->peek() : peek = m_file.peek();
181 if( parsed_event_header && peek=='E' ) break;
182 }
183
184 // Check if all particles in last vertex were parsed
185 /** @bug HepMC2 files produced with Pythia8 are known to have wrong
186 information about number of particles in vertex. Hence '<' sign */
187 if( is_parsing_successful && current_vertex_particles_parsed < current_vertex_particles_count ) {
188 HEPMC3_ERROR( "ReaderAsciiHepMC2: not all particles parsed" )
189 is_parsing_successful = false;
190 }
191 // Check if all vertices were parsed
192 else if( is_parsing_successful && m_vertex_cache.size() != vertices_count ) {
193 HEPMC3_ERROR( "ReaderAsciiHepMC2: not all vertices parsed" )
194 is_parsing_successful = false;
195 }
196
197 if( !is_parsing_successful ) {
198 HEPMC3_ERROR( "ReaderAsciiHepMC2: event parsing failed. Returning empty event" )
199 HEPMC3_DEBUG( 1, "Parsing failed at line:" << std::endl << buf )
200 evt.clear();
201 m_isstream ? m_stream->clear(std::ios::badbit) : m_file.clear(std::ios::badbit);
202 return 0;
203 }
204
205 // Restore production vertex pointers
206 for(unsigned int i=0; i<m_particle_cache.size(); ++i) {
207 if( !m_end_vertex_barcodes[i] ) continue;
208
209 for(unsigned int j=0; j<m_vertex_cache.size(); ++j) {
210 if( m_vertex_barcodes[j] == m_end_vertex_barcodes[i] ) {
211 m_vertex_cache[j]->add_particle_in(m_particle_cache[i]);
212 break;
213 }
214 }
215 }
216
217 // Remove vertices with no incoming particles or no outgoing particles
218 for(unsigned int i=0; i<m_vertex_cache.size(); ++i) {
219 if( m_vertex_cache[i]->particles_in().size() == 0 ) {
220 HEPMC3_DEBUG( 30, "ReaderAsciiHepMC2::read_event - found a vertex without incoming particles: "<<m_vertex_cache[i]->id() );
221//Sometimes the root vertex has no incoming particles. Here we try to save the event.
222 std::vector<GenParticlePtr> beams;
223 for (auto p: m_vertex_cache[i]->particles_out()) if (p->status()==4 && !(p->end_vertex())) beams.push_back(p);
224 for (auto p: beams)
225 {
226 m_vertex_cache[i]->add_particle_in(p);
227 m_vertex_cache[i]->remove_particle_out(p);
228 HEPMC3_DEBUG( 30, "ReaderAsciiHepMC2::read_event - moved particle with status=4 from the outgoing to the incoming particles of vertex: "<<m_vertex_cache[i]->id() );
229 }
230 if (beams.size()==0) m_vertex_cache[i] = nullptr;
231 }
232 else if( m_vertex_cache[i]->particles_out().size() == 0 ) {
233 HEPMC3_DEBUG( 30, "ReaderAsciiHepMC2::read_event - found a vertex without outgouing particles: "<<m_vertex_cache[i]->id() );
234 m_vertex_cache[i] = nullptr;
235 }
236 }
237
238 // Reserve memory for the event
239 evt.reserve( m_particle_cache.size(), m_vertex_cache.size() );
240
241 // Add whole event tree in topological order
242 evt.add_tree( m_particle_cache );
243
244 for(unsigned int i=0; i<m_particle_cache.size(); ++i) {
245 if(m_particle_cache_ghost[i]->attribute_names().size())
246 {
247 std::shared_ptr<DoubleAttribute> phi = m_particle_cache_ghost[i]->attribute<DoubleAttribute>("phi");
248 if (phi) m_particle_cache[i]->add_attribute("phi",phi);
249 std::shared_ptr<DoubleAttribute> theta = m_particle_cache_ghost[i]->attribute<DoubleAttribute>("theta");
250 if (theta) m_particle_cache[i]->add_attribute("theta",theta);
251 if (m_options.find("particle_flows_are_separated")!=m_options.end())
252 {
253 std::shared_ptr<IntAttribute> flow1 = m_particle_cache_ghost[i]->attribute<IntAttribute>("flow1");
254 if (flow1) m_particle_cache[i]->add_attribute("flow1",flow1);
255 std::shared_ptr<IntAttribute> flow2 = m_particle_cache_ghost[i]->attribute<IntAttribute>("flow2");
256 if (flow2) m_particle_cache[i]->add_attribute("flow2",flow2);
257 std::shared_ptr<IntAttribute> flow3 = m_particle_cache_ghost[i]->attribute<IntAttribute>("flow3");
258 if (flow3) m_particle_cache[i]->add_attribute("flow3",flow3);
259 }
260 else
261 {
262 std::shared_ptr<VectorIntAttribute> flows = m_particle_cache_ghost[i]->attribute<VectorIntAttribute>("flows");
263 if (flows) m_particle_cache[i]->add_attribute("flows",flows);
264 }
265 }
266 }
267
268 for(unsigned int i=0; i<m_vertex_cache.size(); ++i)
269 if(m_vertex_cache_ghost[i]->attribute_names().size())
270 {
271 for (size_t ii=0; ii<max_weights_size; ii++)
272 {
273 std::shared_ptr<DoubleAttribute> rs=m_vertex_cache_ghost[i]->attribute<DoubleAttribute>("weight"+std::to_string((long long unsigned int)ii));
274 if (!rs) break;
275 m_vertex_cache[i]->add_attribute("weight"+std::to_string((long long unsigned int)ii),rs);
276 }
277 }
278 std::shared_ptr<IntAttribute> signal_process_vertex_barcode=evt.attribute<IntAttribute>("signal_process_vertex");
279 if (signal_process_vertex_barcode) {
280 int signal_process_vertex_barcode_value=signal_process_vertex_barcode->value();
281 for(unsigned int i=0; i<m_vertex_cache.size(); ++i)
282 {
283 if (i>=m_vertex_barcodes.size()) continue;//this should not happen!
284 if (signal_process_vertex_barcode_value!=m_vertex_barcodes.at(i)) continue;
285 std::shared_ptr<IntAttribute> signal_process_vertex = std::make_shared<IntAttribute>(m_vertex_cache.at(i)->id());
286 evt.add_attribute("signal_process_vertex",signal_process_vertex);
287 break;
288 }
289 }
290 m_particle_cache_ghost.clear();
291 m_vertex_cache_ghost.clear();
292 m_event_ghost->clear();
293 return 1;
294}
295
296int ReaderAsciiHepMC2::parse_event_information(GenEvent &evt, const char *buf) {
297 const char *cursor = buf;
298 int event_no = 0;
299 int vertices_count = 0;
300 int random_states_size = 0;
301 int weights_size = 0;
302 std::vector<long> random_states(0);
303 std::vector<double> weights(0);
304
305 // event number
306 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
307 event_no = atoi(cursor);
308 evt.set_event_number(event_no);
309
310 //mpi
311 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
312 std::shared_ptr<IntAttribute> mpi = std::make_shared<IntAttribute>(atoi(cursor));
313 evt.add_attribute("mpi",mpi);
314
315 //event scale
316 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
317 std::shared_ptr<DoubleAttribute> event_scale = std::make_shared<DoubleAttribute>(atof(cursor));
318 evt.add_attribute("event_scale",event_scale);
319
320 //alpha_qcd
321 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
322 std::shared_ptr<DoubleAttribute> alphaQCD = std::make_shared<DoubleAttribute>(atof(cursor));
323 evt.add_attribute("alphaQCD",alphaQCD);
324
325 //alpha_qed
326 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
327 std::shared_ptr<DoubleAttribute> alphaQED = std::make_shared<DoubleAttribute>(atof(cursor));
328 evt.add_attribute("alphaQED",alphaQED);
329
330 //signal_process_id
331 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
332 std::shared_ptr<IntAttribute> signal_process_id = std::make_shared<IntAttribute>(atoi(cursor));
333 evt.add_attribute("signal_process_id",signal_process_id);
334
335 //signal_process_vertex
336 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
337 std::shared_ptr<IntAttribute> signal_process_vertex = std::make_shared<IntAttribute>(atoi(cursor));
338 evt.add_attribute("signal_process_vertex",signal_process_vertex);
339
340 // num_vertices
341 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
342 vertices_count = atoi(cursor);
343
344 // SKIPPED: beam 1
345 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
346
347 // SKIPPED: beam 2
348 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
349
350 //random states
351 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
352 random_states_size = atoi(cursor);
353 random_states.resize(random_states_size);
354
355 for ( int i = 0; i < random_states_size; ++i ) {
356 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
357 random_states[i] = atoi(cursor);
358 }
359 if (m_options.find("event_random_states_are_separated")!=m_options.end())
360 {
361 evt.add_attribute("random_states",std::make_shared<VectorLongIntAttribute>(random_states));
362 }
363 else
364 {
365 for ( int i = 0; i < random_states_size; ++i )
366 evt.add_attribute("random_states"+std::to_string((long long unsigned int)i),std::make_shared<IntAttribute>(random_states[i]));
367 }
368 // weights
369 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
370 weights_size = atoi(cursor);
371 weights.resize(weights_size);
372
373 for ( int i = 0; i < weights_size; ++i ) {
374 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
375 weights[i] = atof(cursor);
376 }
377
378 evt.weights() = weights;
379
380 HEPMC3_DEBUG( 10, "ReaderAsciiHepMC2: E: "<<event_no<<" ("<<vertices_count<<"V, "<<weights_size<<"W, "<<random_states_size<<"RS)" )
381
382 return vertices_count;
383}
384
385bool ReaderAsciiHepMC2::parse_units(GenEvent &evt, const char *buf) {
386 const char *cursor = buf;
387
388 // momentum
389 if( !(cursor = strchr(cursor+1,' ')) ) return false;
390 ++cursor;
391 Units::MomentumUnit momentum_unit = Units::momentum_unit(cursor);
392
393 // length
394 if( !(cursor = strchr(cursor+1,' ')) ) return false;
395 ++cursor;
396 Units::LengthUnit length_unit = Units::length_unit(cursor);
397
398 evt.set_units(momentum_unit,length_unit);
399
400 HEPMC3_DEBUG( 10, "ReaderAsciiHepMC2: U: " << Units::name(evt.momentum_unit()) << " " << Units::name(evt.length_unit()) )
401
402 return true;
403}
404
405int ReaderAsciiHepMC2::parse_vertex_information(const char *buf) {
406 GenVertexPtr data = std::make_shared<GenVertex>();
407 GenVertexPtr data_ghost = std::make_shared<GenVertex>();
408 FourVector position;
409 const char *cursor = buf;
410 int barcode = 0;
411 int num_particles_out = 0;
412 int weights_size = 0;
413 std::vector<double> weights(0);
414 // barcode
415 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
416 barcode = atoi(cursor);
417
418 // status
419 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
420 data->set_status( atoi(cursor) );
421
422 // x
423 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
424 position.setX(atof(cursor));
425
426 // y
427 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
428 position.setY(atof(cursor));
429
430 // z
431 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
432 position.setZ(atof(cursor));
433
434 // t
435 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
436 position.setT(atof(cursor));
437 data->set_position( position );
438
439 // SKIPPED: num_orphans_in
440 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
441
442 // num_particles_out
443 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
444 num_particles_out = atoi(cursor);
445
446 // weights
447
448 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
449 weights_size = atoi(cursor);
450 weights.resize(weights_size);
451
452 for ( int i = 0; i < weights_size; ++i ) {
453 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
454 weights[i] = atof(cursor);
455 }
456
457
458
459 // Add original vertex barcode to the cache
460 m_vertex_cache.push_back( data );
461 m_vertex_barcodes.push_back( barcode );
462
463 m_event_ghost->add_vertex(data_ghost);
464 if (m_options.find("vertex_weights_are_separated")!=m_options.end())
465 {
466 for ( int i = 0; i < weights_size; ++i )
467 data_ghost->add_attribute("weight"+std::to_string((long long unsigned int)i),std::make_shared<DoubleAttribute>(weights[i]));
468 }
469 else
470 {
471 data_ghost->add_attribute("weights",std::make_shared<VectorDoubleAttribute>(weights));
472 data_ghost->add_attribute("weights",std::make_shared<VectorDoubleAttribute>(weights));
473 }
474 m_vertex_cache_ghost.push_back( data_ghost );
475
476 HEPMC3_DEBUG( 10, "ReaderAsciiHepMC2: V: "<<-(int)m_vertex_cache.size()<<" (old barcode"<<barcode<<") "<<num_particles_out<<" particles)" )
477
478 return num_particles_out;
479}
480
481int ReaderAsciiHepMC2::parse_particle_information(const char *buf) {
482 GenParticlePtr data = std::make_shared<GenParticle>();
483 GenParticlePtr data_ghost = std::make_shared<GenParticle>();
484 m_event_ghost->add_particle(data_ghost);
485 FourVector momentum;
486 const char *cursor = buf;
487 int end_vtx = 0;
488
489 /// @note barcode is ignored
490 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
491
492 // id
493 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
494 data->set_pid( atoi(cursor) );
495
496 // px
497 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
498 momentum.setPx(atof(cursor));
499
500 // py
501 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
502 momentum.setPy(atof(cursor));
503
504 // pz
505 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
506 momentum.setPz(atof(cursor));
507
508 // pe
509 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
510 momentum.setE(atof(cursor));
511 data->set_momentum(momentum);
512
513 // m
514 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
515 data->set_generated_mass( atof(cursor) );
516
517 // status
518 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
519 data->set_status( atoi(cursor) );
520
521 //theta
522 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
523 std::shared_ptr<DoubleAttribute> theta = std::make_shared<DoubleAttribute>(atof(cursor));
524 if (theta->value()!=0.0) data_ghost->add_attribute("theta",theta);
525
526 //phi
527 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
528 std::shared_ptr<DoubleAttribute> phi = std::make_shared<DoubleAttribute>(atof(cursor));
529 if (phi->value()!=0.0) data_ghost->add_attribute("phi",phi);
530
531 // end_vtx_code
532 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
533 end_vtx = atoi(cursor);
534
535 //flow
536 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
537 int flowsize=atoi(cursor);
538
539 std::map<int,int> flows;
540 for (int i=0; i<flowsize; i++)
541 {
542 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
543 int flowindex=atoi(cursor);
544 if( !(cursor = strchr(cursor+1,' ')) ) return -1;
545 int flowvalue=atoi(cursor);
546 flows[flowindex]=flowvalue;
547 }
548 if (m_options.find("particle_flows_are_separated")==m_options.end())
549 {
550 std::vector<int> vectorflows;
551 for (auto f: flows) vectorflows.push_back(f.second);
552 data_ghost->add_attribute("flows",std::make_shared<VectorIntAttribute>(vectorflows));
553 }
554 else
555 {
556 for (auto f: flows) data_ghost->add_attribute("flow"+std::to_string((long long int)f.first),std::make_shared<IntAttribute>(f.second));
557 }
558// Set prod_vtx link
559 if( end_vtx == m_vertex_barcodes.back() ) {
560 m_vertex_cache.back()->add_particle_in(data);
561 end_vtx = 0;
562 }
563 else {
564 m_vertex_cache.back()->add_particle_out(data);
565 }
566
567 m_particle_cache.push_back( data );
568 m_particle_cache_ghost.push_back( data_ghost );
569 m_end_vertex_barcodes.push_back( end_vtx );
570
571 HEPMC3_DEBUG( 10, "ReaderAsciiHepMC2: P: "<<m_particle_cache.size()<<" ( pid: "<<data->pid()<<") end vertex: "<<end_vtx )
572
573 return 0;
574}
575
576bool ReaderAsciiHepMC2::parse_xs_info(GenEvent &evt, const char *buf) {
577 const char *cursor = buf;
578 std::shared_ptr<GenCrossSection> xs = std::make_shared<GenCrossSection>();
579
580 if( !(cursor = strchr(cursor+1,' ')) ) return false;
581 double xs_val = atof(cursor);
582
583 if( !(cursor = strchr(cursor+1,' ')) ) return false;
584 double xs_err = atof(cursor);
585
586 xs->set_cross_section( xs_val, xs_err);
587 evt.add_attribute("GenCrossSection",xs);
588
589 return true;
590}
591
592bool ReaderAsciiHepMC2::parse_weight_names(const char *buf) {
593 const char *cursor = buf;
594 const char *cursor2 = buf;
595 int w_count = 0;
596 std::vector<std::string> w_names;
597
598 // Ignore weight names if no GenRunInfo object
599 if( !run_info() ) return true;
600
601 if( !(cursor = strchr(cursor+1,' ')) ) return false;
602 w_count = atoi(cursor);
603
604 if( w_count <= 0 ) return false;
605
606 w_names.resize(w_count);
607
608 for( int i=0; i < w_count; ++i ) {
609 // Find pair of '"' characters
610 if( !(cursor = strchr(cursor+1,'"')) ) return false;
611 if( !(cursor2 = strchr(cursor+1,'"')) ) return false;
612
613 // Strip cursor of leading '"' character
614 ++cursor;
615
616 w_names[i].assign(cursor, cursor2-cursor);
617
618 cursor = cursor2;
619 }
620
621 run_info()->set_weight_names(w_names);
622
623 return true;
624}
625
626bool ReaderAsciiHepMC2::parse_heavy_ion(GenEvent &evt, const char *buf) {
627 std::shared_ptr<GenHeavyIon> hi = std::make_shared<GenHeavyIon>();
628 const char *cursor = buf;
629
630 if( !(cursor = strchr(cursor+1,' ')) ) return false;
631 hi->Ncoll_hard = atoi(cursor);
632
633 if( !(cursor = strchr(cursor+1,' ')) ) return false;
634 hi->Npart_proj = atoi(cursor);
635
636 if( !(cursor = strchr(cursor+1,' ')) ) return false;
637 hi->Npart_targ = atoi(cursor);
638
639 if( !(cursor = strchr(cursor+1,' ')) ) return false;
640 hi->Ncoll = atoi(cursor);
641
642 if( !(cursor = strchr(cursor+1,' ')) ) return false;
643 hi->spectator_neutrons = atoi(cursor);
644
645 if( !(cursor = strchr(cursor+1,' ')) ) return false;
646 hi->spectator_protons = atoi(cursor);
647
648 if( !(cursor = strchr(cursor+1,' ')) ) return false;
649 hi->N_Nwounded_collisions = atoi(cursor);
650
651 if( !(cursor = strchr(cursor+1,' ')) ) return false;
652 hi->Nwounded_N_collisions = atoi(cursor);
653
654 if( !(cursor = strchr(cursor+1,' ')) ) return false;
655 hi->Nwounded_Nwounded_collisions = atoi(cursor);
656
657 if( !(cursor = strchr(cursor+1,' ')) ) return false;
658 hi->impact_parameter = atof(cursor);
659
660 if( !(cursor = strchr(cursor+1,' ')) ) return false;
661 hi->event_plane_angle = atof(cursor);
662
663 if( !(cursor = strchr(cursor+1,' ')) ) return false;
664 hi->eccentricity = atof(cursor);
665
666 if( !(cursor = strchr(cursor+1,' ')) ) return false;
667 hi->sigma_inel_NN = atof(cursor);
668
669 // Not in HepMC2:
670 hi->centrality = 0.0;
671
672 evt.add_attribute("GenHeavyIon",hi);
673
674 return true;
675}
676
677bool ReaderAsciiHepMC2::parse_pdf_info(GenEvent &evt, const char *buf) {
678 std::shared_ptr<GenPdfInfo> pi = std::make_shared<GenPdfInfo>();
679 const char *cursor = buf;
680
681 if( !(cursor = strchr(cursor+1,' ')) ) return false;
682 pi->parton_id[0] = atoi(cursor);
683
684 if( !(cursor = strchr(cursor+1,' ')) ) return false;
685 pi->parton_id[1] = atoi(cursor);
686
687 if( !(cursor = strchr(cursor+1,' ')) ) return false;
688 pi->x[0] = atof(cursor);
689
690 if( !(cursor = strchr(cursor+1,' ')) ) return false;
691 pi->x[1] = atof(cursor);
692
693 if( !(cursor = strchr(cursor+1,' ')) ) return false;
694 pi->scale = atof(cursor);
695
696 if( !(cursor = strchr(cursor+1,' ')) ) return false;
697 pi->xf[0] = atof(cursor);
698
699 if( !(cursor = strchr(cursor+1,' ')) ) return false;
700 pi->xf[1] = atof(cursor);
701
702 //For compatibility with original HepMC2
703 bool pdfids=true;
704 if( !(cursor = strchr(cursor+1,' ')) ) pdfids=false;
705 if(pdfids) pi->pdf_id[0] = atoi(cursor);
706 else pi->pdf_id[0] =0;
707
708 if(pdfids) if( !(cursor = strchr(cursor+1,' ')) ) pdfids=false;
709 if(pdfids) pi->pdf_id[1] = atoi(cursor);
710 else pi->pdf_id[1] =0;
711
712 evt.add_attribute("GenPdfInfo",pi);
713
714 return true;
715}
716bool ReaderAsciiHepMC2::failed() { return m_isstream ? (bool)m_stream->rdstate() :(bool)m_file.rdstate(); }
717
718void ReaderAsciiHepMC2::close() {
719 if (m_event_ghost) { m_event_ghost->clear(); delete m_event_ghost; m_event_ghost=nullptr;}
720 if( !m_file.is_open() ) return;
721 m_file.close();
722}
723
724} // namespace HepMC3
Note: See TracBrowser for help on using the repository browser.