From 8ce1a0073b25caf9b375552e50f863932eb19bd7 Mon Sep 17 00:00:00 2001
From: Markus Frank <Markus.Frank@cern.ch>
Date: Tue, 13 Feb 2018 21:38:06 +0100
Subject: [PATCH] Fix problem in DDG4 MC truth handling
---
DDCore/include/XML/DocumentHandler.h | 2 +
DDCore/include/XML/XML.h | 1 -
DDCore/include/XML/XMLElements.h | 53 ++++++++++++++++++-
DDCore/src/XML/DocumentHandler.cpp | 22 ++++++--
DDCore/src/XML/XMLElements.cpp | 16 +++++-
DDCore/src/XML/tinyxmlparser_inl.h | 9 ++--
DDG4/include/DDG4/Geant4Particle.h | 1 +
DDG4/src/Geant4Particle.cpp | 68 +++++++++++++++++++++++--
DDG4/src/Geant4ParticleHandler.cpp | 76 ++++++++++++++++++++++------
9 files changed, 214 insertions(+), 34 deletions(-)
diff --git a/DDCore/include/XML/DocumentHandler.h b/DDCore/include/XML/DocumentHandler.h
index ecc5a209c..bb52b1c7a 100644
--- a/DDCore/include/XML/DocumentHandler.h
+++ b/DDCore/include/XML/DocumentHandler.h
@@ -63,6 +63,8 @@ namespace dd4hep {
/// Write xml document to output file (stdout if file name empty)
virtual int output(Document doc, const std::string& fname) const;
+ /// Set minimum print level
+ static int setMinimumPrintLevel(int level);
/// System ID of a given XML entity
static std::string system_path(Handle_t base);
/// System ID of a new XML entity in the same directory as base
diff --git a/DDCore/include/XML/XML.h b/DDCore/include/XML/XML.h
index 64cba6293..d81692608 100644
--- a/DDCore/include/XML/XML.h
+++ b/DDCore/include/XML/XML.h
@@ -33,7 +33,6 @@ typedef dd4hep::xml::DetElement xml_det_t;
typedef dd4hep::xml::Component xml_comp_t;
typedef dd4hep::xml::ChildValue xml_val_t;
typedef dd4hep::xml::Document xml_doc_t;
-typedef dd4hep::xml::Document xml_doc_t;
typedef dd4hep::xml::DocumentHolder xml_doc_holder_t;
typedef dd4hep::xml::DocumentHandler xml_handler_t;
diff --git a/DDCore/include/XML/XMLElements.h b/DDCore/include/XML/XMLElements.h
index d62cfe03f..c88a987dc 100644
--- a/DDCore/include/XML/XMLElements.h
+++ b/DDCore/include/XML/XMLElements.h
@@ -132,8 +132,6 @@ namespace dd4hep {
void _toDictionary(const XmlChar* name, float value);
/// Helper function to populate the evaluator dictionary \ingroup DD4HEP_XML
void _toDictionary(const XmlChar* name, double value);
- /// Helper function to lookup environment from the expression evaluator
- std::string getEnviron(const std::string& env);
/// Conversion function from raw unicode string to bool \ingroup DD4HEP_XML
bool _toBool(const XmlChar* value);
@@ -146,6 +144,11 @@ namespace dd4hep {
/// Conversion function from raw unicode string to double \ingroup DD4HEP_XML
double _toDouble(const XmlChar* value);
+ /// Helper function to lookup environment from the expression evaluator
+ std::string getEnviron(const std::string& env);
+ /// Enable/disable environment resolution when parsing strings
+ bool enableEnvironResolution(bool new_value);
+
/// Helper class to encapsulate a unicode string.
/**
* Simple conversion from ascii strings to unicode strings.
@@ -448,6 +451,8 @@ namespace dd4hep {
void setAttrs(Handle_t e) const;
/// Access typed attribute value by it's unicode name
template <class T> T attr(const XmlChar* name) const;
+ /// Access typed attribute value by it's unicode name. If not existing returns default value
+ template <class T> T attr(const XmlChar* name, T default_value) const;
/// Generic attribute setter with unicode value
Attribute setAttr(const XmlChar* t, const XmlChar* v) const;
@@ -472,6 +477,11 @@ namespace dd4hep {
template <class T> T attr(const char* name) const {
return this->attr<T>(Strng_t(name));
}
+ /// Access typed attribute value by it's name
+ template <class T> T attr(const char* name, const T& default_value) const {
+ Strng_t tag(name);
+ return this->hasAttr(tag) ? this->attr<T>(tag) : default_value;
+ }
/// Generic attribute setter with text value
Attribute setAttr(const XmlChar* t, const char* v) const;
#endif
@@ -537,6 +547,37 @@ namespace dd4hep {
template <> INLINE std::string Handle_t::attr<std::string>(const XmlChar* tag_value) const {
return _toString(attr_value(tag_value));
}
+
+ template <> INLINE bool Handle_t::attr<bool>(const XmlChar* tag_value, bool default_value) const {
+ Attribute a = attr_nothrow(tag_value);
+ return a ? _toBool(attr_value(a)) : default_value;
+ }
+
+ template <> INLINE int Handle_t::attr<int>(const XmlChar* tag_value, int default_value) const {
+ Attribute a = attr_nothrow(tag_value);
+ return a ? _toInt(attr_value(a)) : default_value;
+ }
+
+ template <> INLINE long Handle_t::attr<long>(const XmlChar* tag_value, long default_value) const {
+ Attribute a = attr_nothrow(tag_value);
+ return a ? _toLong(attr_value(a)) : default_value;
+ }
+
+ template <> INLINE float Handle_t::attr<float>(const XmlChar* tag_value, float default_value) const {
+ Attribute a = attr_nothrow(tag_value);
+ return a ? _toFloat(attr_value(a)) : default_value;
+ }
+
+ template <> INLINE double Handle_t::attr<double>(const XmlChar* tag_value, double default_value) const {
+ Attribute a = attr_nothrow(tag_value);
+ return a ? _toDouble(attr_value(a)) : default_value;
+ }
+
+ template <> INLINE std::string Handle_t::attr<std::string>(const XmlChar* tag_value, std::string default_value) const {
+ Attribute a = attr_nothrow(tag_value);
+ return a ? _toString(attr_value(a)) : default_value;
+ }
+
#if 0
template<> INLINE bool Handle_t::attr<bool>(const Attribute tag_value) const
{ return _toBool(attr_value(tag_value));}
@@ -778,11 +819,19 @@ namespace dd4hep {
template <class T> T attr(const XmlChar* tag_value) const {
return m_element.attr<T>(tag_value);
}
+ /// Access attribute with implicit return type conversion
+ template <class T> T attr(const XmlChar* tag_value, T default_value) const {
+ return m_element.attr<T>(tag_value, default_value);
+ }
#ifndef __TIXML__
/// Access typed attribute value by it's name
template <class T> T attr(const char* name) const {
return this->attr<T>(Strng_t(name));
}
+ /// Access typed attribute value by it's name
+ template <class T> T attr(const char* name, T default_value) const {
+ return this->attr<T>(Strng_t(name), default_value);
+ }
#endif
/// Access attribute name (throws exception if not present)
const XmlChar* attr_name(const Attribute a) const {
diff --git a/DDCore/src/XML/DocumentHandler.cpp b/DDCore/src/XML/DocumentHandler.cpp
index 3cd2264ae..3d42e7a95 100644
--- a/DDCore/src/XML/DocumentHandler.cpp
+++ b/DDCore/src/XML/DocumentHandler.cpp
@@ -40,6 +40,7 @@ namespace {
}
return fn;
}
+ int s_minPrintLevel = INFO;
}
#ifndef __TIXML__
@@ -177,8 +178,10 @@ namespace dd4hep {
string baseURI(_toString(id->getBaseURI()));
string schema(_toString(id->getSchemaLocation()));
string ns(_toString(id->getNameSpace()));
- printout(INFO,"XercesC","+++ Resolved URI: sysID:%s uri:%s ns:%s schema:%s",
- systemID.c_str(), baseURI.c_str(), ns.c_str(), schema.c_str());
+ if ( s_minPrintLevel <= INFO ) {
+ printout(INFO,"XercesC","+++ Resolved URI: sysID:%s uri:%s ns:%s schema:%s",
+ systemID.c_str(), baseURI.c_str(), ns.c_str(), schema.c_str());
+ }
#endif
return new MemBufInputSource(input,buf.length(),systemID.c_str(),true);
}
@@ -481,7 +484,7 @@ Document DocumentHandler::load(const std::string& fname, UriReader* reader) cons
printout(WARNING,"DocumentHandler","+++ Loading document URI: %s %s",
fname.c_str(),"[URI Resolution is not supported by TiXML]");
}
- else {
+ else if ( s_minPrintLevel <= INFO ) {
printout(INFO,"DocumentHandler","+++ Loading document URI: %s [Resolved:'%s']",
fname.c_str(),clean.c_str());
}
@@ -506,8 +509,10 @@ Document DocumentHandler::load(const std::string& fname, UriReader* reader) cons
printout(ERROR,"DocumentHandler","+++ Exception (TinyXML): parse(path):%s",e.what());
}
if ( result ) {
- printout(INFO,"DocumentHandler","+++ Document %s succesfully parsed with TinyXML .....",
- fname.c_str());
+ if ( s_minPrintLevel <= INFO ) {
+ printout(INFO,"DocumentHandler","+++ Document %s succesfully parsed with TinyXML .....",
+ fname.c_str());
+ }
return (XmlDocument*)doc;
}
delete doc;
@@ -609,6 +614,13 @@ DocumentHandler::DocumentHandler() {}
/// Default destructor of a document handler using TiXml
DocumentHandler::~DocumentHandler() {}
+/// Set minimum print level
+int DocumentHandler::setMinimumPrintLevel(int level) {
+ int tmp = s_minPrintLevel;
+ s_minPrintLevel = level;
+ return tmp;
+}
+
/// Default comment string
std::string DocumentHandler::defaultComment() {
const char comment[] = "\n"
diff --git a/DDCore/src/XML/XMLElements.cpp b/DDCore/src/XML/XMLElements.cpp
index 821a4a430..5a356494a 100644
--- a/DDCore/src/XML/XMLElements.cpp
+++ b/DDCore/src/XML/XMLElements.cpp
@@ -33,10 +33,15 @@ namespace dd4hep {
}
// Static storage
namespace {
+ bool s_resolve_environment = true;
+
XmlTools::Evaluator& eval(dd4hep::evaluator());
string _checkEnviron(const string& env) {
- string r = getEnviron(env);
- return r.empty() ? env : r;
+ if ( s_resolve_environment ) {
+ string r = getEnviron(env);
+ return r.empty() ? env : r;
+ }
+ return env;
}
}
@@ -424,6 +429,13 @@ string dd4hep::xml::getEnviron(const string& env) {
}
}
+/// Enable/disable environment resolution when parsing strings
+bool dd4hep::xml::enableEnvironResolution(bool new_value) {
+ bool tmp = s_resolve_environment;
+ s_resolve_environment = new_value;
+ return tmp;
+}
+
template <typename B>
static inline string i_add(const string& a, B b) {
string r = a;
diff --git a/DDCore/src/XML/tinyxmlparser_inl.h b/DDCore/src/XML/tinyxmlparser_inl.h
index 9a0ade46b..436d44cb1 100644
--- a/DDCore/src/XML/tinyxmlparser_inl.h
+++ b/DDCore/src/XML/tinyxmlparser_inl.h
@@ -119,20 +119,23 @@ void TiXmlBase::ConvertUTF32ToUTF8( unsigned long input, char* output, int* leng
--output;
*output = (char)((input | BYTE_MARK) & BYTE_MASK);
input >>= 6;
- [[fallthrough]];
+ // fall through [[fallthrough]];
case 3:
--output;
*output = (char)((input | BYTE_MARK) & BYTE_MASK);
input >>= 6;
- [[fallthrough]];
+ // fall through [[fallthrough]];
case 2:
--output;
*output = (char)((input | BYTE_MARK) & BYTE_MASK);
input >>= 6;
- [[fallthrough]];
+ // fall through [[fallthrough]];
case 1:
--output;
*output = (char)(input | FIRST_BYTE_MARK[*length]);
+ // fall through [[fallthrough]];
+ default:
+ break;
}
}
diff --git a/DDG4/include/DDG4/Geant4Particle.h b/DDG4/include/DDG4/Geant4Particle.h
index 9cce168b6..10ec3dd3c 100644
--- a/DDG4/include/DDG4/Geant4Particle.h
+++ b/DDG4/include/DDG4/Geant4Particle.h
@@ -226,6 +226,7 @@ namespace dd4hep {
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 dumpWithMomentumAndVertex(int level, const std::string& src, const char* tag) const;
+ static void header4(int level, const std::string& src, const char* tag);
void dump4(int level, const std::string& src, const char* tag) const;
/// Handlers
diff --git a/DDG4/src/Geant4Particle.cpp b/DDG4/src/Geant4Particle.cpp
index 89aa18c13..61aa060e3 100644
--- a/DDG4/src/Geant4Particle.cpp
+++ b/DDG4/src/Geant4Particle.cpp
@@ -24,6 +24,7 @@
#include "G4ChargedGeantino.hh"
#include "G4Geantino.hh"
+#include <sstream>
#include <iostream>
using namespace dd4hep;
@@ -304,19 +305,75 @@ void Geant4ParticleHandle::dumpWithMomentumAndVertex(int level, const std::strin
text);
}
+void Geant4ParticleHandle::header4(int level, const std::string& src, const char* tag) {
+ printout((dd4hep::PrintLevel)level,src,
+ "+++ %s %10s/%-7s %12s/%-10s %6s/%-6s %4s %4s "
+ "%-4s %-3s %-3s %-10s "
+ "%-5s %-6s %-3s %-3s %-20s %s",
+ tag,"ID", "G4-ID", "Part-Name","PDG", "Parent","G4-ID", "#Par","#Dau",
+ "Prim","Sec",">E","Energy",
+ "EMPTY","STAB","DEC","DOC",
+ "Process", "Processing Flags");
+}
+
void Geant4ParticleHandle::dump4(int level, const std::string& src, const char* tag) const {
Geant4ParticleHandle p(*this);
- char equiv[32];
+ //char equiv[32];
PropertyMask mask(p->reason);
PropertyMask status(p->status);
std::string proc_name = p.processName();
std::string proc_type = p.processTypeName();
+ std::string proc = '['+proc_name+(p->process ? "/" : "")+proc_type+']';
int parent_id = p->parents.empty() ? -1 : *(p->parents.begin());
- equiv[0] = 0;
- if ( p->parents.end() == p->parents.find(p->g4Parent) ) {
- ::snprintf(equiv,sizeof(equiv),"/%d",p->g4Parent);
- }
+ //equiv[0] = 0;
+ //if ( p->parents.end() == p->parents.find(p->g4Parent) ) {
+ // ::snprintf(equiv,sizeof(equiv),"/%d",p->g4Parent);
+ //}
+ std::stringstream str;
+ str << "Parents: ";
+ for( const auto i : p->parents )
+ str << i << " ";
+ str << " Daughters: ";
+ for( const auto i : p->daughters )
+ str << i << " ";
+ printout((dd4hep::PrintLevel)level,src,
+ "+++ %s ID:%7d/%-7d %12s/%-10d %6d/%-6d %4d %4d %-4s %-3s %-3s %+.3e "
+ "%-5s %-4s %-3s %-3s %-20s %c%c%c%c -- %c%c%c%c%c%c%c%c%c %s",
+ tag,
+ p->id,p->originalG4ID,
+ p.particleName().c_str(),
+ p->pdgID,
+ parent_id,p->g4Parent,
+ p.numParent(),
+ int(p->daughters.size()),
+ yes_no(mask.isSet(G4PARTICLE_PRIMARY)),
+ yes_no(mask.isSet(G4PARTICLE_HAS_SECONDARIES)),
+ yes_no(mask.isSet(G4PARTICLE_ABOVE_ENERGY_THRESHOLD)),
+ p.energy(),
+ //
+ yes_no(mask.isSet(G4PARTICLE_CREATED_CALORIMETER_HIT)),
+ yes_no(mask.isSet(G4PARTICLE_CREATED_TRACKER_HIT)),
+ yes_no(mask.isSet(G4PARTICLE_KEEP_PROCESS)),
+ mask.isSet(G4PARTICLE_KEEP_PARENT) ? "YES" : "",
+ proc.c_str(),
+ // 13 flags in total
+ status.isSet(G4PARTICLE_GEN_EMPTY) ? 'E' : '.', // 1
+ status.isSet(G4PARTICLE_GEN_STABLE) ? 'S' : '.',
+ status.isSet(G4PARTICLE_GEN_DECAYED) ? 'D' : '.',
+ status.isSet(G4PARTICLE_GEN_DOCUMENTATION) ? 'd' : '.',
+ status.isSet(G4PARTICLE_GEN_BEAM) ? 'B' : '.', // 5
+ status.isSet(G4PARTICLE_GEN_OTHER) ? 'o' : '.',
+ status.isSet(G4PARTICLE_SIM_CREATED) ? 's' : '.',
+ status.isSet(G4PARTICLE_SIM_BACKSCATTER) ? 'b' : '.',
+ status.isSet(G4PARTICLE_SIM_PARENT_RADIATED) ? 'v' : '.',
+ status.isSet(G4PARTICLE_SIM_DECAY_TRACKER) ? 't' : '.', // 10
+ status.isSet(G4PARTICLE_SIM_DECAY_CALO) ? 'c' : '.',
+ status.isSet(G4PARTICLE_SIM_LEFT_DETECTOR) ? 'l' : '.',
+ status.isSet(G4PARTICLE_SIM_STOPPED) ? 's' : '.',
+ str.str().c_str()
+ );
+#if 0
printout((dd4hep::PrintLevel)level,src,
"+++ %s ID:%7d %12s %6d%-7s %7s %3s %5d %3s %+.3e %-4s %-7s %-3s %-3s %2d [%s%s%s] %c%c%c%c -- %c%c%c%c%c%c%c",
tag,
@@ -351,6 +408,7 @@ void Geant4ParticleHandle::dump4(int level, const std::string& src, const char*
status.isSet(G4PARTICLE_SIM_LEFT_DETECTOR) ? 'l' : '.',
status.isSet(G4PARTICLE_SIM_STOPPED) ? 's' : '.'
);
+#endif
}
/// Default destructor
diff --git a/DDG4/src/Geant4ParticleHandler.cpp b/DDG4/src/Geant4ParticleHandler.cpp
index 6cad82846..f12a745e3 100644
--- a/DDG4/src/Geant4ParticleHandler.cpp
+++ b/DDG4/src/Geant4ParticleHandler.cpp
@@ -306,7 +306,6 @@ void Geant4ParticleHandler::end(const G4Track* track) {
ph->vey = p.y();
ph->vez = p.z();
-
// Set the simulator status bits
PropertyMask simStatus(m_currTrack.status);
@@ -388,8 +387,10 @@ void Geant4ParticleHandler::beginEvent(const G4Event* event) {
/// Debugging: Dump Geant4 particle map
void Geant4ParticleHandler::dumpMap(const char* tag) const {
+ const string& n = name();
+ Geant4ParticleHandle::header4(INFO,n,tag);
for(ParticleMap::const_iterator iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i) {
- Geant4ParticleHandle((*i).second).dump4(INFO,name(),tag);
+ Geant4ParticleHandle((*i).second).dump4(INFO,n,tag);
}
}
@@ -398,14 +399,14 @@ void Geant4ParticleHandler::endEvent(const G4Event* event) {
int count = 0;
int level = outputLevel();
do {
- if ( level <= VERBOSE ) dumpMap("Particle");
+ if ( level <= VERBOSE ) dumpMap("Particle ");
debug("+++ Iteration:%d Tracks:%d Equivalents:%d",++count,m_particleMap.size(),m_equivalentTracks.size());
} while( recombineParents() > 0 );
- if ( level <= VERBOSE ) dumpMap("Recombined");
+ if ( level <= VERBOSE ) dumpMap( "Recombined");
// Rebase the simulated tracks, so that they fit to the generator particles
rebaseSimulatedTracks(0);
- if ( level <= VERBOSE ) dumpMap("Rebased");
+ if ( level <= VERBOSE ) dumpMap( "Rebased ");
// Consistency check....
checkConsistency();
/// Call the user particle handler
@@ -466,8 +467,8 @@ void Geant4ParticleHandler::rebaseSimulatedTracks(int ) {
}
TrackEquivalents::mapped_type equiv = (*ie).second;
if ( ipar != m_particleMap.end() ) {
- equivalents[(*ie).first] = (*ipar).second->id; // requires (1) !
Geant4ParticleHandle p = (*ipar).second;
+ equivalents[(*ie).first] = p->id; // requires (1) to be filled properly!
const G4ParticleDefinition* def = p.definition();
int pdg = int(fabs(def->GetPDGEncoding())+0.1);
if ( pdg != 0 && pdg<36 && !(pdg > 10 && pdg < 17) && pdg != 22 ) {
@@ -486,22 +487,65 @@ void Geant4ParticleHandler::rebaseSimulatedTracks(int ) {
// (3) Compute the particle's parents and daughters.
// Replace the original Geant4 track with the
// equivalent particle still present in the record.
- for(iend=m_particleMap.end(), i=m_particleMap.begin(); i!=iend; ++i) {
+ // Note:
+ // We rely here on the ordering of the particles accoding to their
+ // Processing by Geant4 to establish mother daughter relationships.
+ // == > use finalParticles map and NOT m_particleMap.
+ int equiv_id = -1;
+ for(iend=finalParticles.end(), i=finalParticles.begin(); i!=iend; ++i) {
Particle* p = (*i).second;
if ( p->g4Parent > 0 ) {
- int equiv_id = equivalents[p->g4Parent];
- if ( (ipar=finalParticles.find(equiv_id)) != finalParticles.end() ) {
- Particle* q = (*ipar).second;
- q->daughters.insert(p->id);
- p->parents.insert(q->id);
+ TrackEquivalents::iterator iequ = equivalents.find(p->g4Parent);
+ if ( iequ != equivalents.end() ) {
+ equiv_id = (*iequ).second;//equivalents[p->g4Parent];
+ if ( (ipar=finalParticles.find(equiv_id)) != finalParticles.end() ) {
+ Particle* q = (*ipar).second;
+ bool prim = (p->reason&G4PARTICLE_PRIMARY) == G4PARTICLE_PRIMARY;
+ // We assume that the mother daughter relationship
+ // is filled by the event readers!
+ if ( !prim ) {
+ p->parents.insert(q->id);
+ }
+ if ( !p->parents.empty() ) {
+ int parent_id = (*p->parents.begin());
+ if ( parent_id == q->id )
+ q->daughters.insert(p->id);
+ else if ( !prim )
+ error("+++ Inconsistency in equivalent record! Parent: %d Daughter:%d",q->id, p->id);
+ }
+ else {
+ error("+++ Inconsistency in parent relashionship: %d NO parent!", p->id);
+ }
+ continue;
+ }
}
- else {
- error("+++ Inconsistency in particle record: Geant4 parent %d "
- "of particle %d (equiv:%d) not in record!",
- p->g4Parent,p->id,equiv_id);
+ error("+++ Inconsistency in particle record: Geant4 parent %d "
+ "of particle %d not in record of final particles!",
+ p->g4Parent,p->id);
+ }
+ }
+#if 0
+ for(iend=finalParticles.end(), i=finalParticles.begin(); i!=iend; ++i) {
+ Particle* p = (*i).second;
+ if ( p->g4Parent > 0 ) {
+ int parent_id = (*p->parents.begin());
+ if ( (ipar=finalParticles.find(parent_id)) != finalParticles.end() ) {
+ Particle* q = (*ipar).second;
+ // Generator particles have a proper history.
+ // We only deal with particles, which are not of MC origin.
+ //p->parents.insert(q->id);
+ if ( parent_id == q->id )
+ q->daughters.insert(p->id);
+ else
+ error("+++ Inconsistency in equivalent record! Parent: %d Daughter:%d",q->id, p->id);
+ continue;
}
+ error("+++ Inconsistency in particle record: Geant4 parent %d "
+ "of particle %d not in record of final particles!",
+ p->g4Parent,p->id);
}
}
+#endif
m_equivalentTracks = equivalents;
m_particleMap = finalParticles;
}
--
GitLab