Newer
Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
//====================================================================
// AIDA Detector description implementation for LCD
//--------------------------------------------------------------------
//
//====================================================================
// Framework include files
#include "DD4hep/IoStreams.h"
#include "DDG4/Geant4InputAction.h"
// C/C++ include files
#include <fstream>
#include <boost/iostreams/stream.hpp>
/// Namespace for the AIDA detector description toolkit
namespace DD4hep {
/// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
namespace Simulation {
namespace HepMC { class EventStream; }
/// Class to populate Geant4 primaries from StdHep files.
/**
* Class to populate Geant4 primary particles and vertices from a
* file in HEPEvt format (ASCII)
*
* \author P.Kostka (main author)
* \author M.Frank (code reshuffeling into new DDG4 scheme)
* \version 1.0
* \ingroup DD4HEP_SIMULATION
*/
class Geant4EventReaderHepMC : public Geant4EventReader {
typedef boost::iostreams::stream<dd4hep_file_source<int> > in_stream;
//typedef boost::iostreams::stream<dd4hep_file_source<TFile*> > in_stream;
typedef HepMC::EventStream EventStream;
protected:
in_stream m_input;
EventStream* m_events;
public:
/// Initializing constructor
Geant4EventReaderHepMC(const std::string& nam);
/// Default destructor
virtual ~Geant4EventReaderHepMC();
/// Read an event and fill a vector of MCParticles.
virtual EventReaderStatus readParticles(int event_number, std::vector<Particle*>& particles);
};
} /* End namespace Simulation */
} /* End namespace DD4hep */
// $Id: Geant4Converter.cpp 603 2013-06-13 21:15:14Z markus.frank $
//====================================================================
// AIDA Detector description implementation for LCD
//--------------------------------------------------------------------
//
//====================================================================
// #include "DDG4/Geant4EventReaderHepMC"
#include "DDG4/Factories.h"
#include "DD4hep/Printout.h"
#include "DDG4/Geant4Primary.h"
#include "CLHEP/Units/SystemOfUnits.h"
#include "CLHEP/Units/PhysicalConstants.h"
#include <cerrno>
#include <algorithm>
using namespace std;
using namespace CLHEP;
using namespace DD4hep::Simulation;
typedef DD4hep::ReferenceBitMask<int> PropertyMask;
// Factory entry
DECLARE_GEANT4_EVENT_READER(Geant4EventReaderHepMC)
/// Namespace for the AIDA detector description toolkit
namespace DD4hep {
/// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit
namespace Simulation {
namespace HepMC {
class EventHeader {
public:
int id;
int num_vertices;
int bp1, bp2;
int signal_process_id;
int signal_process_vertex;
float scale;
float alpha_qcd;
float alpha_qed;
vector<float> weights;
vector<long> random;
EventHeader() {}
};
/// The known_io enum is used to track which type of input is being read
enum known_io { gen=1, ascii, extascii, ascii_pdt, extascii_pdt };
class EventStream {
public:
typedef std::map<int,Geant4Vertex*> Vertices;
typedef std::map<int,Geant4Particle*> Particles;
// io information
string key;
double mom_unit, pos_unit;
int io_type;
float xsection, xsection_err;
EventHeader header;
Vertices m_vertices;
Particles m_particles;
EventStream(istream& in) : instream(in), io_type(0)
{ use_default_units(); }
/// Check if data stream is in proper state and has data
bool ok() const;
Geant4Vertex* vertex(int i);
Particles& particles() { return m_particles; }
Vertices& vertices() { return m_vertices; }
void set_io(int typ, const string& k)
{ io_type = typ; key = k; }
void use_default_units()
bool read();
void clear();
};
char get_input(istream& is, istringstream& iline);
int read_weight_names(EventStream &, istringstream& iline);
int read_particle(EventStream &info, istringstream& iline, Geant4Particle * p);
int read_vertex(EventStream &info, istream& is, istringstream & iline);
int read_event_header(EventStream &info, istringstream & input, EventHeader& header);
int read_cross_section(EventStream &info, istringstream & input);
int read_units(EventStream &info, istringstream & input);
int read_heavy_ion(EventStream &, istringstream & input);
int read_pdf(EventStream &, istringstream & input);
Geant4Vertex* vertex(EventStream& info, int i);
void fix_particles(EventStream &info);
}
}
}
/// Initializing constructor
Geant4EventReaderHepMC::Geant4EventReaderHepMC(const string& nam)
: Geant4EventReader(nam), m_input(), m_events(0)
{
// Now open the input file:
m_input.open(nam.c_str(),BOOST_IOS::in|BOOST_IOS::binary);
if ( m_input < 0 ) {
string err = "+++ Geant4EventReaderHepMC: Failed to open input stream:"+nam+
" Error:"+string(strerror(errno));
throw runtime_error(err);
}
m_events = new HepMC::EventStream(m_input);
}
/// Default destructor
Geant4EventReaderHepMC::~Geant4EventReaderHepMC() {
delete m_events;
m_events = 0;
m_input.close();
}
/// Read an event and fill a vector of MCParticles.
Geant4EventReaderHepMC::EventReaderStatus
Geant4EventReaderHepMC::readParticles(int /* ev_id */, Particles& output) {
return EVENT_READER_IO_ERROR;
}
EventStream::Particles& p = m_events->particles();
output.reserve(p.size());
transform(p.begin(),p.end(),back_inserter(output),reference2nd(p));
m_events->clear();
for(Particles::const_iterator k=output.begin(); k != output.end(); ++k) {
Geant4ParticleHandle p(*k);
printout(INFO,m_name,
"+++ %s ID:%3d status:%08X typ:%9d Mom:(%+.2e,%+.2e,%+.2e)[MeV] "
"time: %+.2e [ns] #Dau:%3d #Par:%1d",
"",p->id,p->status,p->pdgID,
p->psx/MeV,p->psy/MeV,p->psz/MeV,p->time/ns,
p->daughters.size(),
p->parents.size());
//output.push_back(p);
}
return EVENT_READER_OK;
}
return EVENT_READER_IO_ERROR;
}
void HepMC::fix_particles(EventStream& info) {
EventStream::Particles& parts = info.particles();
EventStream::Vertices& verts = info.vertices();
EventStream::Particles::iterator i;
std::set<int>::const_iterator id, ip;
Geant4ParticleHandle p((*i).second);
int end_vtx_id = p->secondaries;
p->secondaries = 0;
if ( v ) {
p->vex = v->x;
p->vey = v->y;
p->vez = v->z;
v->in.insert(p->id);
for(id=v->out.begin(); id!=v->out.end();++id)
p->daughters.insert(*id);
}
}
EventStream::Vertices::iterator j;
Geant4Vertex* v = (*j).second;
for (id=v->in.begin(); id!=v->in.end();++id) {
for (ip=v->out.begin(); ip!=v->out.end();++ip) {
EventStream::Particles::iterator ipp = parts.find(*ip);
(*ipp).second->parents.insert(*id);
}
}
}
}
Geant4Vertex* HepMC::vertex(EventStream& info, int i) {
EventStream::Vertices::iterator it=info.vertices().find(i);
return (it==info.vertices().end()) ? 0 : (*it).second;
}
char HepMC::get_input(istream& is, istringstream& iline) {
char value = is.peek();
if ( !is ) { // make sure the stream is valid
cerr << "StreamHelpers: setting badbit." << endl;
is.clear(ios::badbit);
return -1;
}
string line, firstc;
getline(is,line);
if ( !is ) { // make sure the stream is valid
cerr << "StreamHelpers: setting badbit." << endl;
is.clear(ios::badbit);
return -1;
}
iline.str(line);
if ( line.length()>0 && line[1] == ' ' ) iline >> firstc;
return iline ? value : -1;
}
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
string line;
while ( is ) {
char val = is.peek();
if( val == 'E' ) { // next event
return 1;
}
getline(is,line);
if ( !is.good() ) return 0;
}
return 0;
}
int HepMC::read_weight_names(EventStream&, istringstream&) {
#if 0
int HepMC::read_weight_names(EventStream& info, istringstream& iline)
size_t name_size = 0;
iline >> name_size;
info.weights.names.clear();
info.weights.weights.clear();
string name;
WeightContainer namedWeight;
string::size_type i1 = line.find("\""), i2, len = line.size();
for(size_t ii = 0; ii < name_size; ++ii) {
// weight names may contain blanks
if(i1 >= len) {
cout << "debug: attempting to read past the end of the named weight line " << endl;
cout << "debug: We should never get here" << endl;
cout << "debug: Looking for the end of this event" << endl;
}
i2 = line.find("\"",i1+1);
name = line.substr(i1+1,i2-i1-1);
if ( namedWeight.names.find(name) == namedWeight.names.end() )
namedWeight.names[name] = namedWeight.names.size();
namedWeight.weights[ii] = info.weights.weights[ii];
i1 = line.find("\"",i2+1);
}
info.weights.weights = namedWeight.weights;
info.weights.names = namedWeight.names;
#endif
return 1;
}
int HepMC::read_particle(EventStream &info, istringstream& input, Geant4Particle * p) {
float ene = 0., theta = 0., phi = 0;
int size = 0;
// check that the input stream is still OK after reading item
input >> p->id >> p->pdgID >> p->psx >> p->psy >> p->psz >> ene;
p->psx /= info.mom_unit;
p->psy /= info.mom_unit;
p->psz /= info.mom_unit;
ene /= info.mom_unit;
if ( !input )
else if ( info.io_type != ascii ) {
input >> p->mass;
p->mass /= info.mom_unit;
}
else {
p->mass = sqrt(fabs(ene*ene - p->psx*p->psx + p->psy*p->psy + p->psz*p->psz));
// Reuse here the secondaries to store the end-vertex ID
input >> p->status >> theta >> phi >> p->secondaries >> size;
if(!input)
return 0;
// read flow patterns if any exist
for (int i = 0; i < size; ++i ) {
input >> p->colorFlow[0] >> p->colorFlow[1];
if(!input) return 0;
}
return 1;
}
int HepMC::read_vertex(EventStream &info, istream& is, istringstream & input) {
int id=0, dummy=0, num_orphans_in=0, num_particles_out=0, weights_size=0;
vector<float> weights;
Geant4Vertex* v = new Geant4Vertex();
Geant4Particle* p;
input >> id >> dummy >> v->x >> v->y >> v->z >> v->time
>> num_orphans_in >> num_particles_out >> weights_size;
v->x /= info.pos_unit;
v->y /= info.pos_unit;
v->z /= info.pos_unit;
weights.resize(weights_size);
for (int i1 = 0; i1 < weights_size; ++i1) {
//cout << "Add Vertex:" << id << endl;
for(char value = is.peek(); value=='P'; value=is.peek()) {
value = get_input(is,input);
if( !input || value < 0 )
read_particle(info, input,p = new Geant4Particle());
if(!input) {
cerr << "Failed to read particle!" << endl;
delete p;
return 0;
}
p->id = info.particles().size();
info.particles().insert(make_pair(p->id,p));
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
p->pex = p->psx;
p->pey = p->psy;
p->pez = p->psz;
if ( --num_orphans_in >= 0 ) {
p->vex = v->x;
p->vey = v->y;
p->vez = v->z;
v->in.insert(p->id);
//cout << "Add INGOING Particle:" << p->id << endl;
}
else if ( num_particles_out >= 0 ) {
v->out.insert(p->id);
p->vsx = v->x;
p->vsy = v->y;
p->vsz = v->z;
//cout << "Add OUTGOING Particle:" << p->id << endl;
}
else {
delete p;
throw runtime_error("Invalid number of particles....");
}
}
return 1;
}
int HepMC::read_event_header(EventStream &info, istringstream & input, EventHeader& header) {
// read values into temp variables, then fill GenEvent
int random_states_size = 0, nmpi = -1;
input >> header.id;
if( info.io_type == gen || info.io_type == extascii ) {
input >> nmpi;
if( input.fail() ) return 0;
//MSF set_mpi( nmpi );
}
input >> header.scale;
input >> header.alpha_qcd;
input >> header.alpha_qed;
input >> header.signal_process_id;
input >> header.signal_process_vertex;
input >> header.num_vertices;
if( info.io_type == gen || info.io_type == extascii )
input >> header.bp1 >> header.bp2;
input >> random_states_size;
cout << input.str() << endl;
input.clear();
if( input.fail() ) return 0;
header.random.resize(random_states_size);
for(int i = 0; i < random_states_size; ++i )
input >> header.random[i];
size_t weights_size = 0;
input >> weights_size;
if( input.fail() ) return 0;
vector<float> wgt(weights_size);
for(size_t ii = 0; ii < weights_size; ++ii )
input >> wgt[ii];
if( input.fail() ) return 0;
// weight names will be added later if they exist
if( weights_size > 0 ) header.weights = wgt;
return 1;
}
int HepMC::read_cross_section(EventStream &info, istringstream & input) {
input >> info.xsection >> info.xsection_err;
return input.fail() ? 0 : 1;
}
int HepMC::read_units(EventStream &info, istringstream & input) {
if( info.io_type == gen ) {
string mom, pos;
input >> mom >> pos;
if ( !input.fail() ) {
if ( mom == "KEV" ) info.mom_unit = keV;
else if ( mom == "MEV" ) info.mom_unit = MeV;
else if ( mom == "GEV" ) info.mom_unit = GeV;
else if ( mom == "TEV" ) info.mom_unit = TeV;
if ( pos == "MM" ) info.pos_unit = mm;
else if ( pos == "CM" ) info.pos_unit = cm;
else if ( pos == "M" ) info.pos_unit = m;
}
}
return input.fail() ? 0 : 1;
}
int HepMC::read_heavy_ion(EventStream &, istringstream & input) {
// read values into temp variables, then create a new HeavyIon object
int nh =0, np =0, nt =0, nc =0,
neut = 0, prot = 0, nw =0, nwn =0, nwnw =0;
float impact = 0., plane = 0., xcen = 0., inel = 0.;
input >> nh >> np >> nt >> nc >> neut >> prot >> nw >> nwn >> nwnw;
input >> impact >> plane >> xcen >> inel;
cerr << "Reading heavy ion, but igoring data!" << endl;
/*
ion->set_Ncoll_hard(nh);
ion->set_Npart_proj(np);
ion->set_Npart_targ(nt);
ion->set_Ncoll(nc);
ion->set_spectator_neutrons(neut);
ion->set_spectator_protons(prot);
ion->set_N_Nwounded_collisions(nw);
ion->set_Nwounded_N_collisions(nwn);
ion->set_Nwounded_Nwounded_collisions(nwnw);
ion->set_impact_parameter(impact);
ion->set_event_plane_angle(plane);
ion->set_eccentricity(xcen);
ion->set_sigma_inel_NN(inel);
*/
return input.fail() ? 0 : 1;
}
int HepMC::read_pdf(EventStream &, istringstream & input) {
// read values into temp variables, then create a new PdfInfo object
int id1 =0, id2 =0, pdf_id1=0, pdf_id2=0;
double x1 = 0., x2 = 0., scale = 0., pdf1 = 0., pdf2 = 0.;
input >> id1 ;
if ( input.fail() )
return 0;
// check now for empty PdfInfo line
if( id1 == 0 )
return 0;
// continue reading
input >> id2 >> x1 >> x2 >> scale >> pdf1 >> pdf2;
if ( input.fail() )
return 0;
// check to see if we are at the end of the line
if( !input.eof() ) {
input >> pdf_id1 >> pdf_id2;
}
cerr << "Reading pdf, but igoring data!" << endl;
/*
pdf->set_id1( id1 );
pdf->set_id2( id2 );
pdf->set_pdf_id1( pdf_id1 );
pdf->set_pdf_id2( pdf_id2 );
pdf->set_x1( x1 );
pdf->set_x2( x2 );
pdf->set_scalePDF( scale );
pdf->set_pdf1( pdf1 );
pdf->set_pdf2( pdf2 );
*/
return input.fail() ? 0 : 1;
}
/// Check if data stream is in proper state and has data
bool HepMC::EventStream::ok() const {
// make sure the stream is good
if ( instream.eof() || instream.fail() ) {
instream.clear(ios::badbit);
return false;
}
return true;
}
void HepMC::EventStream::clear() {
releaseObjects(m_vertices)();
releaseObjects(m_particles)();
}
bool HepMC::EventStream::read() {
EventStream& info = *this;
releaseObjects(vertices())();
releaseObjects(particles())();
while( instream.good() ) {
char value = instream.peek();
istringstream input_line;
if ( value == 'E' && event_read ) break;
else if ( instream.eof() ) return false;
else if ( value=='#' || ::isspace(value) ) {
get_input(instream,input_line);
continue;
}
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
value = get_input(instream,input_line);
// On failure switch to end
if( !input_line || value < 0 )
goto Skip;
switch( value ) {
case 'H': {
int iotype = 0;
string key;
input_line >> key;
// search for event listing key before first event only.
key = key.substr(0,key.find('\r'));
if ( key == "H" && (this->io_type == gen || this->io_type == extascii) ) {
read_heavy_ion(info, input_line);
break;
}
else if( key == "HepMC::IO_GenEvent-START_EVENT_LISTING" )
this->set_io(gen,key);
else if( key == "HepMC::IO_Ascii-START_EVENT_LISTING" )
this->set_io(ascii,key);
else if( key == "HepMC::IO_ExtendedAscii-START_EVENT_LISTING" )
this->set_io(extascii,key);
else if( key == "HepMC::IO_Ascii-START_PARTICLE_DATA" )
this->set_io(ascii_pdt,key);
else if( key == "HepMC::IO_ExtendedAscii-START_PARTICLE_DATA" )
this->set_io(extascii_pdt,key);
else if( key == "HepMC::IO_GenEvent-END_EVENT_LISTING" )
iotype = gen;
else if( key == "HepMC::IO_Ascii-END_EVENT_LISTING" )
iotype = ascii;
else if( key == "HepMC::IO_ExtendedAscii-END_EVENT_LISTING" )
iotype = extascii;
else if( key == "HepMC::IO_Ascii-END_PARTICLE_DATA" )
iotype = ascii_pdt;
else if( key == "HepMC::IO_ExtendedAscii-END_PARTICLE_DATA" )
iotype = extascii_pdt;
if( iotype != 0 && this->io_type != iotype ) {
cerr << "GenEvent::find_end_key: iotype keys have changed. "
<< "MALFORMED INPUT" << endl;
instream.clear(ios::badbit);
return false;
}
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
continue;
}
case 'E': // deal with the event line
if ( !read_event_header(info, input_line, this->header) )
goto Skip;
event_read = true;
continue;
case 'N': // get weight names
if ( !read_weight_names(info, input_line) )
goto Skip;
continue;
case 'U': // get unit information if it exists
if ( !read_units(info, input_line) )
goto Skip;
continue;
case 'C': // we have a GenCrossSection line
if ( !read_cross_section(info, input_line) )
goto Skip;
continue;
case 'V': // Read vertex with particles
if ( !read_vertex(info, instream, input_line) )
goto Skip;
continue;
case 'F': // Read PDF
if ( !read_pdf(info, input_line) )
goto Skip;
continue;
case 'P': // we should not find this line
cerr << "streaming input: found unexpected Particle line." << endl;
continue;
default: // ignore everything else
continue;
}
continue;
Skip:
printout(WARNING,"HepMC::EventStream","+++ Skip event with ID: %d",this->header.id);
releaseObjects(vertices())();
releaseObjects(particles())();
event_read = false;
if ( instream.eof() ) return false;
}
releaseObjects(vertices())();
return true;
}