From 8fee842fc750c8bd1a2f1a2538dec49d43b556d0 Mon Sep 17 00:00:00 2001
From: Markus Frank <markus.frank@cern.ch>
Date: Mon, 27 Oct 2014 13:20:03 +0000
Subject: [PATCH] Add first version of HepMC file reader. Move IoStreams back
to DDG4
---
DDCore/include/DD4hep/Primitives.h | 3 +
DDG4/examples/CLICSidSimu.py | 7 ++-
DDG4/include/DDG4/Geant4Particle.h | 3 +-
.../DD4hep => DDG4/include/DDG4}/IoStreams.h | 0
DDG4/plugins/Geant4EventReaderHepEvt.cpp | 9 +--
DDG4/plugins/Geant4EventReaderHepMC.cpp | 62 ++++++++++++++-----
DDG4/src/Geant4InputAction.cpp | 9 +--
DDG4/src/Geant4InputHandling.cpp | 8 +--
DDG4/src/Geant4Particle.cpp | 25 +++++++-
DDG4/src/Geant4ParticleGun.cpp | 2 +-
DDG4/src/Geant4ParticleHandler.cpp | 5 +-
{DDCore => DDG4}/src/IoStreams.cpp | 2 +-
12 files changed, 97 insertions(+), 38 deletions(-)
rename {DDCore/include/DD4hep => DDG4/include/DDG4}/IoStreams.h (100%)
rename {DDCore => DDG4}/src/IoStreams.cpp (99%)
diff --git a/DDCore/include/DD4hep/Primitives.h b/DDCore/include/DD4hep/Primitives.h
index 719de23fa..ecebfac5a 100644
--- a/DDCore/include/DD4hep/Primitives.h
+++ b/DDCore/include/DD4hep/Primitives.h
@@ -300,6 +300,9 @@ namespace DD4hep {
void clear(const T& m) {
mask &= ~m;
}
+ void clear() {
+ mask = 0;
+ }
bool isSet(const T& m) const {
return (mask&m) == m;
}
diff --git a/DDG4/examples/CLICSidSimu.py b/DDG4/examples/CLICSidSimu.py
index 425fffb96..4642ca2be 100644
--- a/DDG4/examples/CLICSidSimu.py
+++ b/DDG4/examples/CLICSidSimu.py
@@ -110,12 +110,13 @@ def run():
#gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/smuonLR_PointK_3TeV_BS_noBkg_run0001.stdhep"
#gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/bbbb_3TeV.stdhep"
#gen.Input = "LCIOFileReader|/home/frankm/SW/data/mcparticles_pi-_5GeV.slcio"
- gen.Input = "LCIOFileReader|/home/frankm/SW/data/mcparticles_mu-_5GeV.slcio"
+ #gen.Input = "LCIOFileReader|/home/frankm/SW/data/mcparticles_mu-_5GeV.slcio"
#gen.Input = "LCIOFileReader|/home/frankm/SW/data/bbbb_3TeV.slcio"
#gen.Input = "LCIOStdHepReader|/home/frankm/SW/data/FCC-eh.stdhep"
- #gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/data.hepmc.txt"
+ gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/data.hepmc.txt"
+ gen.Input = "Geant4EventReaderHepMC|/home/frankm/SW/data/sherpa-2.1.1_zjets.hepmc2g"
gen.OutputLevel = 4 # generator_output_level
- gen.MomentumScale = 0.1
+ gen.MomentumScale = 1.0
gen.Mask = 1
gen.enableUI()
kernel.generatorAction().adopt(gen)
diff --git a/DDG4/include/DDG4/Geant4Particle.h b/DDG4/include/DDG4/Geant4Particle.h
index 6e31d970a..ed3614f86 100644
--- a/DDG4/include/DDG4/Geant4Particle.h
+++ b/DDG4/include/DDG4/Geant4Particle.h
@@ -196,7 +196,8 @@ namespace DD4hep {
/// Output type 2:+++ "tag" 20 G4: 7 def:0xde4eaa8 [gamma , gamma] reason: 20 E:+3.304035e+01 in record:YES \#Par: 1/18 \#Dau: 0
void dump2(int level, const std::string& src, const char* tag, int g4id, bool inrec) const;
/// Output type 3:+++ "tag" ID: 0 e- status:00000014 type: 11 Vertex:(+0.00e+00,+0.00e+00,+0.00e+00) [mm] time: +0.00e+00 [ns] \#Par: 0 \#Dau: 4
- void dump3(int level, const std::string& src, const char* tag) const;
+ void dumpWithVertex(int level, const std::string& src, const char* tag) const;
+ void dumpWithMomentum(int level, const std::string& src, const char* tag) const;
void dump4(int level, const std::string& src, const char* tag) const;
/// Handlers
diff --git a/DDCore/include/DD4hep/IoStreams.h b/DDG4/include/DDG4/IoStreams.h
similarity index 100%
rename from DDCore/include/DD4hep/IoStreams.h
rename to DDG4/include/DDG4/IoStreams.h
diff --git a/DDG4/plugins/Geant4EventReaderHepEvt.cpp b/DDG4/plugins/Geant4EventReaderHepEvt.cpp
index ff4c3e053..cb87b858d 100644
--- a/DDG4/plugins/Geant4EventReaderHepEvt.cpp
+++ b/DDG4/plugins/Geant4EventReaderHepEvt.cpp
@@ -192,10 +192,11 @@ Geant4EventReaderHepEvt::readParticles(int /* event_number */, vector<Particle*>
// Generator status
// Simulator status 0 until simulator acts on it
p->status = 0;
- if ( ISTHEP == 0 ) status.set(G4PARTICLE_GEN_EMPTY);
- if ( ISTHEP == 1 ) status.set(G4PARTICLE_GEN_STABLE);
- if ( ISTHEP == 2 ) status.set(G4PARTICLE_GEN_DECAYED);
- if ( ISTHEP == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
+ if ( ISTHEP == 0 ) status.set(G4PARTICLE_GEN_EMPTY);
+ else if ( ISTHEP == 1 ) status.set(G4PARTICLE_GEN_STABLE);
+ else if ( ISTHEP == 2 ) status.set(G4PARTICLE_GEN_DECAYED);
+ else if ( ISTHEP == 3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
+ else status.set(G4PARTICLE_GEN_DOCUMENTATION);
//
// Creation time (note the units [1/c_light])
// (No information in HEPEvt files)
diff --git a/DDG4/plugins/Geant4EventReaderHepMC.cpp b/DDG4/plugins/Geant4EventReaderHepMC.cpp
index 62d95d535..d34df09b4 100644
--- a/DDG4/plugins/Geant4EventReaderHepMC.cpp
+++ b/DDG4/plugins/Geant4EventReaderHepMC.cpp
@@ -218,13 +218,31 @@ void HepMC::fix_particles(EventStream& info) {
EventStream::Vertices::iterator j;
for(j=verts.begin(); j != verts.end(); ++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);
+ for (ip=v->out.begin(); ip!=v->out.end();++ip) {
+ EventStream::Particles::iterator ipp = parts.find(*ip);
+ Geant4Particle* p = (*ipp).second;
+ for (id=v->in.begin(); id!=v->in.end();++id) {
+ p->parents.insert(*id);
}
}
}
+ /// Particles originating from the beam (=no parents) must be
+ /// be stripped off their parents and the status set to G4PARTICLE_GEN_DECAYED!
+ vector<Geant4Particle*> beam;
+ for(i=parts.begin(); i != parts.end(); ++i) {
+ Geant4ParticleHandle p((*i).second);
+ if ( p->parents.size() == 0 ) {
+ for(id=p->daughters.begin(); id!=p->daughters.end();++id) {
+ Geant4Particle *pp = parts[*id];
+ beam.push_back(pp);
+ }
+ }
+ }
+ for(vector<Geant4Particle*>::iterator i=beam.begin(); i!=beam.end();++i) {
+ cout << "Clear parents of " << (*i)->id << endl;
+ (*i)->parents.clear();
+ (*i)->status = G4PARTICLE_GEN_DECAYED;
+ }
}
Geant4Vertex* HepMC::vertex(EventStream& info, int i) {
@@ -299,28 +317,41 @@ int HepMC::read_weight_names(EventStream&, istringstream&) {
int HepMC::read_particle(EventStream &info, istringstream& input, Geant4Particle * p) {
float ene = 0., theta = 0., phi = 0;
- int size = 0;
+ int size = 0, stat=0;
+ PropertyMask status(p->status);
// 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;
+ p->id = info.particles().size();
+ p->psx *= info.mom_unit;
+ p->psy *= info.mom_unit;
+ p->psz *= info.mom_unit;
+ ene *= info.mom_unit;
if ( !input )
return 0;
else if ( info.io_type != ascii ) {
input >> p->mass;
- p->mass /= info.mom_unit;
+ 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;
+ input >> stat >> theta >> phi >> p->secondaries >> size;
if(!input)
return 0;
+ //
+ // Generator status
+ // Simulator status 0 until simulator acts on it
+ status.clear();
+ if ( stat == 0 ) status.set(G4PARTICLE_GEN_EMPTY);
+ else if ( stat == 0x1 ) status.set(G4PARTICLE_GEN_STABLE);
+ else if ( stat == 0x2 ) status.set(G4PARTICLE_GEN_DECAYED);
+ else if ( stat == 0x3 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
+ else if ( stat == 0x4 ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
+ else if ( stat == 0xB ) status.set(G4PARTICLE_GEN_DOCUMENTATION);
+
// read flow patterns if any exist
for (int i = 0; i < size; ++i ) {
input >> p->colorFlow[0] >> p->colorFlow[1];
@@ -339,9 +370,9 @@ int HepMC::read_vertex(EventStream &info, istream& is, istringstream & input)
>> num_orphans_in >> num_particles_out >> weights_size;
if(!input)
return 0;
- v->x /= info.pos_unit;
- v->y /= info.pos_unit;
- v->z /= info.pos_unit;
+ 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) {
input >> weights[i1];
@@ -362,16 +393,15 @@ int HepMC::read_vertex(EventStream &info, istream& is, istringstream & input)
delete p;
return 0;
}
- p->id = info.particles().size();
info.particles().insert(make_pair(p->id,p));
p->pex = p->psx;
p->pey = p->psy;
p->pez = p->psz;
if ( --num_orphans_in >= 0 ) {
+ v->in.insert(p->id);
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 ) {
diff --git a/DDG4/src/Geant4InputAction.cpp b/DDG4/src/Geant4InputAction.cpp
index e5c095576..bfbf58ef9 100644
--- a/DDG4/src/Geant4InputAction.cpp
+++ b/DDG4/src/Geant4InputAction.cpp
@@ -159,9 +159,10 @@ void Geant4InputAction::operator()(G4Event* event) {
Geant4ParticleHandle p(primaries[i]);
const double mom_scale = m_momScale;
PropertyMask status(p->status);
- p->psx = mom_scale*p->psx;
- p->psy = mom_scale*p->psy;
- p->psz = mom_scale*p->psz;
+ p->psx = mom_scale*p->psx;
+ p->psy = mom_scale*p->psy;
+ p->psz = mom_scale*p->psz;
+
if ( p->parents.size() == 0 ) {
if ( status.isSet(G4PARTICLE_GEN_EMPTY) || status.isSet(G4PARTICLE_GEN_DOCUMENTATION) )
vtx->in.insert(p->id); // Beam particles and primary quarks etc.
@@ -169,6 +170,6 @@ void Geant4InputAction::operator()(G4Event* event) {
vtx->out.insert(p->id); // Stuff, to be given to Geant4 together with daughters
}
inter->particles.insert(make_pair(p->id,p));
- p.dump3(outputLevel()-1,name(),"+->");
+ p.dumpWithVertex(outputLevel()-1,name(),"+->");
}
}
diff --git a/DDG4/src/Geant4InputHandling.cpp b/DDG4/src/Geant4InputHandling.cpp
index baf1a422f..2e09c3a0d 100644
--- a/DDG4/src/Geant4InputHandling.cpp
+++ b/DDG4/src/Geant4InputHandling.cpp
@@ -301,6 +301,7 @@ int DD4hep::Simulation::generatePrimaries(const Geant4Action* caller,
Interaction::VertexMap& vm = interaction->vertices;
map<int,G4PrimaryParticle*> prim;
set<int> visited;
+ char text[64];
Geant4PrimaryInteraction::VertexMap::iterator ivfnd, iv, ivend;
for(Interaction::VertexMap::const_iterator iend=vm.end(),i=vm.begin(); i!=iend; ++i) {
@@ -324,11 +325,8 @@ int DD4hep::Simulation::generatePrimaries(const Geant4Action* caller,
PropertyMask reason(r->reason);
reason.set(G4PARTICLE_PRIMARY);
v4->SetPrimary(p4);
- caller->printM1("+++ +-> G4Primary[%3d] ID:%3d type:%9d/%-12s "
- "Momentum:(%+.2e,%+.2e,%+.2e) [GeV] time:%+.2e [ns] #Par:%3d #Dau:%3d",
- num_part,r->id,r->pdgID,r.particleName().c_str(),
- r->psx/GeV,r->psy/GeV,r->psz/GeV,r->time/ns,
- int(r->parents.size()),int(r->daughters.size()));
+ ::snprintf(text,sizeof(text),"+-> G4Primary[%3d]",num_part);
+ r.dumpWithMomentum(caller->outputLevel()-1,caller->name(),text);
++num_part;
}
}
diff --git a/DDG4/src/Geant4Particle.cpp b/DDG4/src/Geant4Particle.cpp
index fd1feb806..656b95710 100644
--- a/DDG4/src/Geant4Particle.cpp
+++ b/DDG4/src/Geant4Particle.cpp
@@ -257,7 +257,7 @@ void Geant4ParticleHandle::dump2(int level, const std::string& src, const char*
}
/// Output type 3:+++ <tag> ID: 0 e- status:00000014 type: 11 Vertex:(+0.00e+00,+0.00e+00,+0.00e+00) [mm] time: +0.00e+00 [ns] \#Par: 0 \#Dau: 4
-void Geant4ParticleHandle::dump3(int level, const std::string& src, const char* tag) const {
+void Geant4ParticleHandle::dumpWithVertex(int level, const std::string& src, const char* tag) const {
char text[256];
Geant4ParticleHandle p(*this);
text[0]=0;
@@ -278,6 +278,29 @@ void Geant4ParticleHandle::dump3(int level, const std::string& src, const char*
text);
}
+
+/// Output type 3:+++ <tag> ID: 0 e- status:00000014 type: 11 Vertex:(+0.00e+00,+0.00e+00,+0.00e+00) [mm] time: +0.00e+00 [ns] \#Par: 0 \#Dau: 4
+void Geant4ParticleHandle::dumpWithMomentum(int level, const std::string& src, const char* tag) const {
+ char text[256];
+ Geant4ParticleHandle p(*this);
+ text[0]=0;
+ if ( p->parents.size() == 1 )
+ ::snprintf(text,sizeof(text),"/%d",*(p->parents.begin()));
+ else if ( p->parents.size() > 1 ) {
+ text[0]='/';text[1]=0;
+ for(std::set<int>::const_iterator i=p->parents.begin(); i!=p->parents.end(); ++i)
+ ::snprintf(text+strlen(text),sizeof(text)-strlen(text),"%d ",*i);
+ }
+ printout((DD4hep::PrintLevel)level,src,
+ "+++ %s ID:%3d %-12s status:%08X typ:%9d Mom:(%+.2e,%+.2e,%+.2e)[mm] "
+ "time: %+.2e [ns] #Dau:%3d #Par:%1d%-6s",
+ tag,p->id,p.particleName().c_str(),p->status,p->pdgID,
+ p->psx/MeV,p->psy/MeV,p->psz/MeV,p->time/ns,
+ int(p->daughters.size()),
+ int(p->parents.size()),
+ text);
+}
+
void Geant4ParticleHandle::dump4(int level, const std::string& src, const char* tag) const {
Geant4ParticleHandle p(*this);
char equiv[32];
diff --git a/DDG4/src/Geant4ParticleGun.cpp b/DDG4/src/Geant4ParticleGun.cpp
index 216ec2d49..d2ec80adf 100644
--- a/DDG4/src/Geant4ParticleGun.cpp
+++ b/DDG4/src/Geant4ParticleGun.cpp
@@ -119,7 +119,7 @@ void Geant4ParticleGun::operator()(G4Event* event) {
status.set(G4PARTICLE_GEN_STABLE);
vtx->out.insert(p->id);
inter->particles.insert(make_pair(p->id,p));
- p.dump3(outputLevel()-1,name(),"+->");
+ p.dumpWithVertex(outputLevel()-1,name(),"+->");
}
++m_shotNo;
diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp
index c2b3155d7..c39f4b905 100644
--- a/DDG4/src/Geant4ParticleHandler.cpp
+++ b/DDG4/src/Geant4ParticleHandler.cpp
@@ -343,7 +343,7 @@ void Geant4ParticleHandler::end(const G4Track* track) {
if ( ip != m_particleMap.end() )
(*ip).second->reason |= track_reason;
else
- ph.dump3(outputLevel()+3,name(),"FATAL: No real particle parent present");
+ ph.dumpWithVertex(outputLevel()+3,name(),"FATAL: No real particle parent present");
}
}
@@ -594,7 +594,7 @@ void Geant4ParticleHandler::checkConsistency() const {
/// First check the consistency of the particle map itself
for(ParticleMap::const_iterator j, i=m_particleMap.begin(); i!=m_particleMap.end(); ++i) {
- Particle* p = (*i).second;
+ Geant4ParticleHandle p((*i).second);
PropertyMask mask(p->reason);
PropertyMask status(p->status);
set<int>& daughters = p->daughters;
@@ -621,6 +621,7 @@ void Geant4ParticleHandler::checkConsistency() const {
char parent_list[1024];
parent_list[0] = 0;
++num_errors;
+ p.dumpWithMomentum(ERROR,name(),"INCONSISTENCY");
for(set<int>::const_iterator ip=p->parents.begin(); ip!=p->parents.end();++ip)
::snprintf(parent_list+strlen(parent_list),sizeof(parent_list)-strlen(parent_list),"%d ",*ip);
error("+++ Particle:%d Parent %d (G4id:%d) In record:%s In parent list:%s [%s]",
diff --git a/DDCore/src/IoStreams.cpp b/DDG4/src/IoStreams.cpp
similarity index 99%
rename from DDCore/src/IoStreams.cpp
rename to DDG4/src/IoStreams.cpp
index fcb90d192..c18cff50e 100644
--- a/DDCore/src/IoStreams.cpp
+++ b/DDG4/src/IoStreams.cpp
@@ -8,7 +8,7 @@
//====================================================================
// Framework includes
-#include "DD4hep/IoStreams.h"
+#include "DDG4/IoStreams.h"
#include <sys/types.h>
#include <sys/stat.h>
--
GitLab