From d3947be8251947e0939841f220235b0ddb0f40c1 Mon Sep 17 00:00:00 2001 From: Markus Frank <markus.frank@cern.ch> Date: Tue, 28 Jun 2016 10:20:16 +0000 Subject: [PATCH] Add multi-segmentations and readouts. See doc/release.notes --- DDCore/include/DD4hep/Alignment.h | 2 +- DDCore/include/DD4hep/Dictionary.h | 2 + DDCore/include/DD4hep/IDDescriptor.h | 10 +- DDCore/include/DD4hep/Readout.h | 10 + .../include/DD4hep/objects/DetectorInterna.h | 4 +- .../include/DD4hep/objects/ObjectsInterna.h | 28 ++ DDCore/include/XML/UnicodeValues.h | 8 + DDCore/include/XML/XMLDimension.h | 18 +- DDCore/src/LCDDImp.cpp | 1 + DDCore/src/Readout.cpp | 58 ++- DDCore/src/Segmentations.cpp | 1 + DDCore/src/VolumeManager.cpp | 180 ++++++++- DDCore/src/XML/XMLDimension.cpp | 5 + DDCore/src/plugins/Compact2Objects.cpp | 163 ++++++-- .../src/CylindricalEndcapCalorimeter_geo.cpp | 2 + DDG4/include/DDG4/Geant4Action.h | 3 + DDG4/include/DDG4/Geant4ReadoutVolumeFilter.h | 61 +++ DDG4/include/DDG4/Geant4SensDetAction.h | 38 +- DDG4/include/DDG4/Geant4SensDetAction.inl | 51 ++- DDG4/plugins/Geant4SDActions.cpp | 18 +- DDG4/plugins/Geant4TrackerWeightedSD.cpp | 2 +- DDG4/python/DDG4.py | 71 +++- DDG4/src/Geant4OutputAction.cpp | 4 +- DDG4/src/Geant4ReadoutVolumeFilter.cpp | 60 +++ DDG4/src/Geant4SensDetAction.cpp | 23 +- .../DDSegmentation/MultiSegmentation.h | 90 +++++ .../include/DDSegmentation/Segmentation.h | 4 +- DDSegmentation/src/MultiSegmentation.cpp | 91 +++++ DDSegmentation/src/Segmentation.cpp | 367 +++++++++--------- cmake/DD4hepConfig.cmake.in | 4 +- doc/release.notes | 16 + examples/CLICSiD/compact/compact.xml | 4 +- examples/ClientTests/CMakeLists.txt | 18 +- .../ClientTests/compact/MultiCollections.xml | 64 +++ .../compact/MultiSegmentCollections.xml | 67 ++++ .../compact/MultiSegmentations.xml | 62 +++ .../ClientTests/scripts/MultiCollections.py | 62 +++ 37 files changed, 1390 insertions(+), 282 deletions(-) create mode 100644 DDG4/include/DDG4/Geant4ReadoutVolumeFilter.h create mode 100644 DDG4/src/Geant4ReadoutVolumeFilter.cpp create mode 100644 DDSegmentation/include/DDSegmentation/MultiSegmentation.h create mode 100644 DDSegmentation/src/MultiSegmentation.cpp create mode 100644 examples/ClientTests/compact/MultiCollections.xml create mode 100644 examples/ClientTests/compact/MultiSegmentCollections.xml create mode 100644 examples/ClientTests/compact/MultiSegmentations.xml create mode 100644 examples/ClientTests/scripts/MultiCollections.py diff --git a/DDCore/include/DD4hep/Alignment.h b/DDCore/include/DD4hep/Alignment.h index 20b48c321..08a59d653 100644 --- a/DDCore/include/DD4hep/Alignment.h +++ b/DDCore/include/DD4hep/Alignment.h @@ -1,4 +1,4 @@ -// $Id: $ +// $Id:$ //========================================================================== // AIDA Detector description implementation for LCD //-------------------------------------------------------------------------- diff --git a/DDCore/include/DD4hep/Dictionary.h b/DDCore/include/DD4hep/Dictionary.h index f9fb5e6bf..7ea791c3f 100644 --- a/DDCore/include/DD4hep/Dictionary.h +++ b/DDCore/include/DD4hep/Dictionary.h @@ -179,9 +179,11 @@ template class DD4hep::Handle<TNamed>; #pragma link C++ class DD4hep::Geometry::Segmentation+; #pragma link C++ class DD4hep::Geometry::SegmentationObject+; #pragma link C++ class DD4hep::Handle<DD4hep::Geometry::SegmentationObject>+; +#pragma link C++ class DD4hep::Geometry::HitCollection+; #pragma link C++ class DD4hep::Geometry::Readout+; #pragma link C++ class DD4hep::Geometry::ReadoutObject+; #pragma link C++ class DD4hep::Handle<DD4hep::Geometry::ReadoutObject>+; +#pragma link C++ class vector<DD4hep::Geometry::HitCollection>+; #pragma link C++ class vector<DD4hep::Geometry::Readout>+; #pragma link C++ class vector<DD4hep::Geometry::IDDescriptor>+; diff --git a/DDCore/include/DD4hep/IDDescriptor.h b/DDCore/include/DD4hep/IDDescriptor.h index 0930f8a78..f98ad9304 100644 --- a/DDCore/include/DD4hep/IDDescriptor.h +++ b/DDCore/include/DD4hep/IDDescriptor.h @@ -44,12 +44,12 @@ namespace DD4hep { */ class IDDescriptor: public Handle<IDDescriptorObject> { public: - typedef IDDescriptorObject Object; - typedef BitFieldValue* Field; - typedef std::vector<std::pair<std::string, Field> > FieldMap; + typedef IDDescriptorObject Object; + typedef BitFieldValue* Field; + typedef std::vector<std::pair<std::string, Field> > FieldMap; typedef std::vector<std::pair<size_t, std::string> > FieldIDs; - typedef std::pair<Field, VolumeID> VolIDField; - typedef std::vector<VolIDField> VolIDFields; + typedef std::pair<Field, VolumeID> VolIDField; + typedef std::vector<VolIDField> VolIDFields; public: /// Default constructor diff --git a/DDCore/include/DD4hep/Readout.h b/DDCore/include/DD4hep/Readout.h index 3998ee6aa..cdeba3c36 100644 --- a/DDCore/include/DD4hep/Readout.h +++ b/DDCore/include/DD4hep/Readout.h @@ -33,9 +33,14 @@ namespace DD4hep { // Forward declarations class DetElement; class ReadoutObject; + class HitCollection; /// Handle to the implementation of the readout structure of a subdetector /** + * If there is no explicit hit collection defined, by default one single + * hit collection is defined by the name of the readout itself. + * If hit collections are defined, ALL must be defined. + * * \author M.Frank * \version 1.0 * \ingroup DD4HEP_GEOMETRY @@ -44,6 +49,7 @@ namespace DD4hep { public: /// Implementation type typedef ReadoutObject Object; + typedef HitCollection Collection; public: /// Default constructor Readout() @@ -71,6 +77,10 @@ namespace DD4hep { m_element = ro.m_element; return *this; } + /// Access explicit names of hit collections if present + std::vector<std::string> collectionNames() const; + /// Access hit collections if present + std::vector<const HitCollection*> collections() const; /// Assign IDDescription to readout structure void setIDDescriptor(const Ref_t& spec) const; /// Access IDDescription structure diff --git a/DDCore/include/DD4hep/objects/DetectorInterna.h b/DDCore/include/DD4hep/objects/DetectorInterna.h index be7c3f3fa..d73644b6f 100644 --- a/DDCore/include/DD4hep/objects/DetectorInterna.h +++ b/DDCore/include/DD4hep/objects/DetectorInterna.h @@ -90,6 +90,8 @@ namespace DD4hep { HAVE_WORLD_TRAFO = 1<<0, HAVE_PARENT_TRAFO = 1<<1, HAVE_REFERENCE_TRAFO = 1<<2, + HAVE_SENSITIVE_DETECTOR = 1<<29, + IS_TOP_LEVEL_DETECTOR = 1<<30, HAVE_OTHER = 1<<31 }; @@ -101,7 +103,7 @@ namespace DD4hep { int id; /// Flag to process hits int combineHits; - /// Flag to encode detector types + /// Flag to encode detector types unsigned int typeFlag; /// Full path to this detector element. May be invalid std::string path; diff --git a/DDCore/include/DD4hep/objects/ObjectsInterna.h b/DDCore/include/DD4hep/objects/ObjectsInterna.h index 5cd8a4583..d75eb9c60 100644 --- a/DDCore/include/DD4hep/objects/ObjectsInterna.h +++ b/DDCore/include/DD4hep/objects/ObjectsInterna.h @@ -135,6 +135,30 @@ namespace DD4hep { virtual ~LimitSetObject(); }; + /// Definition of the HitCollection parameters used by the Readout + /** + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_GEOMETRY + */ + class HitCollection { + public: + /// Hit collection name + std::string name; + /// Discriminator key name from the <id/> string + std::string key; + /// Range values of the key is not empty. + long key_min, key_max; + /// Default constructor + HitCollection() {} + /// Copy constructor + HitCollection(const HitCollection& c); + /// Initializing constructor + HitCollection(const std::string& name, const std::string& key="",long k_min=~0x0, long kmax=~0x0); + /// Assignment operator + HitCollection& operator=(const HitCollection& c); + }; + /// Concrete object implementation of the Readout Handle /** * @@ -144,12 +168,16 @@ namespace DD4hep { */ class ReadoutObject: public NamedObject { public: + typedef HitCollection Collection; + typedef std::vector<HitCollection> Collections; /// Handle to the readout segmentation Segmentation segmentation; /// Handle to the volume Volume readoutWorld; /// Handle to the field descriptor IDDescriptor id; + /// Hit collection container (if defined) + Collections hits; /// Standard constructor ReadoutObject(); /// Default destructor diff --git a/DDCore/include/XML/UnicodeValues.h b/DDCore/include/XML/UnicodeValues.h index 5be6e1769..70a40dcd8 100644 --- a/DDCore/include/XML/UnicodeValues.h +++ b/DDCore/include/XML/UnicodeValues.h @@ -77,6 +77,8 @@ namespace DD4hep { UNICODE (coefficient); UNICODE (coefficients); UNICODE (color); + UNICODE (collections); + UNICODE (collection); UNICODE (combine_hits); UNICODE (combineHits); UNICODE (comment); @@ -173,6 +175,7 @@ namespace DD4hep { UNICODE (half_z); UNICODE (header); UNICODE (height); + UNICODE (hits_collections); UNICODE (hits_collection); UNICODE (hype); @@ -205,6 +208,11 @@ namespace DD4hep { UNICODE (J); UNICODE (k); + UNICODE (key); + UNICODE (key_min); + UNICODE (key_max); + UNICODE (key_val); + UNICODE (key_value); UNICODE (K); UNICODE (l); diff --git a/DDCore/include/XML/XMLDimension.h b/DDCore/include/XML/XMLDimension.h index 7ce5fafb6..cd1d9eac1 100644 --- a/DDCore/include/XML/XMLDimension.h +++ b/DDCore/include/XML/XMLDimension.h @@ -467,14 +467,24 @@ namespace DD4hep { /// Access attribute values: phi_tilt double phi_tilt() const; /// Access attribute values: nphi - int nphi() const; + int nphi() const; /// Access attribute values: rc double rc() const; /// Access attribute values: zstart double zstart() const; /// Access attribute values: nz - int nz() const; + int nz() const; + /// Access attribute values: key + int key() const; + /// Access attribute values: key_min + int key_min() const; + /// Access attribute values: key_max + int key_max() const; + /// Access attribute values: key_val + int key_val() const; + /// Access attribute values: key_value + int key_value() const; /// Access attribute values: start double start() const; @@ -490,9 +500,9 @@ namespace DD4hep { double outer_field() const; /// Access attribute values: visible - bool visible() const; + bool visible() const; /// Access attribute values: show_daughters - bool show_daughters() const; + bool show_daughters() const; /// Access min/max parameters: cut double cut() const; diff --git a/DDCore/src/LCDDImp.cpp b/DDCore/src/LCDDImp.cpp index 6d642bbfe..b116a8c79 100644 --- a/DDCore/src/LCDDImp.cpp +++ b/DDCore/src/LCDDImp.cpp @@ -217,6 +217,7 @@ LCDD& LCDDImp::addDetector(const Ref_t& ref_det) { } } m_detectors.append(ref_det); + det_element->flag |= DetElement::Object::IS_TOP_LEVEL_DETECTOR; Volume volume = det_element.placement()->GetMotherVolume(); if ( volume == m_worldVol ) { printout(DEBUG,"LCDD","Added detector %s to the world instance.",det_element.name()); diff --git a/DDCore/src/Readout.cpp b/DDCore/src/Readout.cpp index 23705c8d3..c1b4875d8 100644 --- a/DDCore/src/Readout.cpp +++ b/DDCore/src/Readout.cpp @@ -25,14 +25,65 @@ using namespace DD4hep; using namespace DD4hep::Geometry; using dd4hep::mm; +/// Copy constructor +HitCollection::HitCollection(const HitCollection& c) + : name(c.name), key(c.key), key_min(c.key_min), key_max(c.key_max) +{ +} + +/// Initializing constructor +HitCollection::HitCollection(const std::string& n, const std::string& k, long k_min, long k_max) + : name(n), key(k), key_min(k_min), key_max(k_max) +{ +} + +/// Assignment operator +HitCollection& HitCollection::operator=(const HitCollection& c) { + if ( this != &c ) { + name = c.name; + key = c.key; + key_min = c.key_min; + key_max = c.key_max; + } + return *this; +} + /// Initializing constructor to create a new object Readout::Readout(const string& nam) { assign(new ReadoutObject(), nam, "readout"); } +/// Access names of hit collections +std::vector<std::string> Readout::collectionNames() const { + std::vector<std::string> colls; + if ( isValid() ) { + Object& ro = object<Object>(); + if ( !ro.hits.empty() ) { + for(Object::Collections::const_iterator i=ro.hits.begin(); i!=ro.hits.end(); ++i) + colls.push_back((*i).name); + } + return colls; + } + throw runtime_error("DD4hep: Readout::collectionsNames: Cannot access object data [Invalid Handle]"); +} + +/// Access hit collections +std::vector<const HitCollection*> Readout::collections() const { + std::vector<const HitCollection*> colls; + if ( isValid() ) { + Object& ro = object<Object>(); + if ( !ro.hits.empty() ) { + for(Object::Collections::const_iterator i=ro.hits.begin(); i!=ro.hits.end(); ++i) + colls.push_back(&(*i)); + } + return colls; + } + throw runtime_error("DD4hep: Readout::collectionsNames: Cannot access object data [Invalid Handle]"); +} + /// Assign IDDescription to readout structure void Readout::setIDDescriptor(const Ref_t& new_descriptor) const { - if (isValid()) { // Remember: segmentation is NOT owned by readout structure! + if ( isValid() ) { // Remember: segmentation is NOT owned by readout structure! if (new_descriptor.isValid()) { // Do NOT delete! data<Object>()->id = new_descriptor; Segmentation seg = data<Object>()->segmentation; @@ -53,10 +104,10 @@ IDDescriptor Readout::idSpec() const { /// Assign segmentation structure to readout void Readout::setSegmentation(const Segmentation& seg) const { - if (isValid()) { + if ( isValid() ) { Object& ro = object<Object>(); Segmentation::Implementation* e = ro.segmentation.ptr(); - if (e) { // Remember: segmentation is owned by readout structure! + if (e) { // Remember: segmentation is owned by readout structure! delete e; // Need to delete the segmentation object } if (seg.isValid()) { @@ -71,3 +122,4 @@ void Readout::setSegmentation(const Segmentation& seg) const { Segmentation Readout::segmentation() const { return object<Object>().segmentation; } + diff --git a/DDCore/src/Segmentations.cpp b/DDCore/src/Segmentations.cpp index 99f8aa77f..f666fd02a 100644 --- a/DDCore/src/Segmentations.cpp +++ b/DDCore/src/Segmentations.cpp @@ -108,6 +108,7 @@ Segmentation::Segmentation(const string& typ, const string& nam) : BaseSegmentation* s = DDSegmentation::SegmentationFactory::instance()->create(typ); if (s != 0) { assign(new Object(s), nam, ""); + if ( !nam.empty() ) s->setName(nam); } else { throw runtime_error("FAILED to create segmentation: " + typ + ". Missing factory method for: " + typ + "!"); } diff --git a/DDCore/src/VolumeManager.cpp b/DDCore/src/VolumeManager.cpp index 8f0b5fa31..69b38d75b 100644 --- a/DDCore/src/VolumeManager.cpp +++ b/DDCore/src/VolumeManager.cpp @@ -33,6 +33,7 @@ namespace { struct Populator { typedef PlacedVolume::VolIDs VolIDs; typedef vector<TGeoNode*> Chain; + typedef pair<VolumeID, VolumeID> Encoding; /// Reference to the LCDD instance LCDD& m_lcdd; /// Reference to the volume manager to be populated @@ -47,14 +48,22 @@ namespace { /// Populate the Volume manager void populate(DetElement e) { const DetElement::Children& c = e.children(); + SensitiveDetector parent_sd; + if ( e->flag&DetElement::Object::HAVE_SENSITIVE_DETECTOR ) { + parent_sd = m_lcdd.sensitiveDetector(e.name()); + } for (DetElement::Children::const_iterator i = c.begin(); i != c.end(); ++i) { DetElement de = (*i).second; PlacedVolume pv = de.placement(); if (pv.isValid()) { VolIDs ids; Chain chain; - SensitiveDetector sd; + SensitiveDetector sd = parent_sd; m_entries.clear(); +#if 0 + Encoding coding; + scanPhysicalVolume(de, de, pv, coding, ids, sd, chain); +#endif scanPhysicalVolume(de, de, pv, ids, sd, chain); continue; } @@ -80,6 +89,96 @@ namespace { return DetElement(); } /// Scan a single physical volume and look for sensitive elements below + size_t scanPhysicalVolume(DetElement& parent, DetElement e, PlacedVolume pv, + Encoding parent_encoding, VolIDs ids, SensitiveDetector& sd, Chain& chain) + { + TGeoNode* node = pv.ptr(); + size_t count = 0; + if (node) { + bool got_readout = false; + Volume vol = pv.volume(); + const VolIDs& pv_ids = pv.volIDs(); + Encoding volume_encoding = parent_encoding; + if ( (parent->flag&DetElement::Object::HAVE_SENSITIVE_DETECTOR) ) { + sd = m_lcdd.sensitiveDetector(parent.name()); + } + else if ( (e->flag&DetElement::Object::HAVE_SENSITIVE_DETECTOR) ) { + sd = m_lcdd.sensitiveDetector(e.name()); + } + else if ( vol.isSensitive() ) { + sd = vol.sensitiveDetector(); + } + chain.push_back(node); + if ( sd.isValid() ) { + Readout ro = sd.readout(); + if ( ro.isValid() ) { + got_readout = true; + if ( ids.empty() ) + volume_encoding = update_encoding(ro.idSpec(), pv_ids, parent_encoding); + else { + ids.VolIDs::Base::insert(ids.end(), pv_ids.begin(), pv_ids.end()); + volume_encoding = update_encoding(ro.idSpec(), ids, parent_encoding); + ids.clear(); + } + add_entry(sd, parent, e, node, volume_encoding, chain); + ++count; + } + else { + printout(WARNING, "VolumeManager", + "%s: Strange constellation volume %s is sensitive, but has no readout! sd:%p", + parent.name(), pv.volume().name(), sd.ptr()); + } + } + if ( !got_readout && !pv_ids.empty() ) { + ids.VolIDs::Base::insert(ids.end(), pv_ids.begin(), pv_ids.end()); + } + for (Int_t idau = 0, ndau = node->GetNdaughters(); idau < ndau; ++idau) { + TGeoNode* daughter = node->GetDaughter(idau); + PlacedVolume placement(daughter); + if ( placement.data() ) { + size_t cnt; + PlacedVolume pv_dau = Ref_t(daughter); + DetElement de_dau = findElt(e, daughter); + if ( de_dau.isValid() ) { + Chain dau_chain; + cnt = scanPhysicalVolume(parent, de_dau, pv_dau, volume_encoding, ids, sd, dau_chain); + } + else { + cnt = scanPhysicalVolume(parent, e, pv_dau, volume_encoding, ids, sd, chain); + } + // There was a sensitive daughter volume, also add the parent entry. + if ( !got_readout && count == 0 && cnt > 0 && sd.isValid() && !pv_ids.empty()) { + add_entry(sd, parent, e, node, volume_encoding, chain); + } + count += cnt; + } + else { + throw runtime_error("Invalid not instrumented placement:"+string(daughter->GetName())+ + " [Internal error -- bad detector constructor]"); + } + } + if ( count == 0 ) { + sd = SensitiveDetector(0); + } + else if ( count > 0 && sd.isValid() ) { + // We recuperate volumes from lower levels by reusing the subdetector + // This only works if there is exactly one sensitive detector per subdetector! + // I hate this, but I could not talk Frank out of this! M.F. + Readout ro = sd.readout(); + if ( ro.isValid() ) { + IDDescriptor iddesc = ro.idSpec(); + Encoding det_encoding = encoding(iddesc,ids); + printout(VERBOSE,"VolumeManager","++++ %-11s SD:%s VolID=%p Mask=%p",e.path().c_str(), + got_readout ? "RECUPERATED" : "REGULAR", sd.name(), + (void*)det_encoding.first, (void*)det_encoding.second); + e.object<DetElement::Object>().volumeID = det_encoding.first; + } + } + chain.pop_back(); + } + return count; + } + /// Scan a single physical volume and look for sensitive elements below size_t scanPhysicalVolume(DetElement& parent, DetElement e, PlacedVolume pv, VolIDs ids, SensitiveDetector& sd, Chain& chain) { @@ -140,7 +239,7 @@ namespace { Readout ro = sd.readout(); if ( ro.isValid() ) { IDDescriptor iddesc = ro.idSpec(); - pair<VolumeID, VolumeID> det_encoding = encoding(iddesc,ids); + Encoding det_encoding = encoding(iddesc,ids); printout(VERBOSE,"VolumeManager","++++ %-11s SD:%s VolID=%p Mask=%p",e.path().c_str(), got_readout ? "RECUPERATED" : "REGULAR", sd.name(), (void*)det_encoding.first, (void*)det_encoding.second); @@ -153,7 +252,20 @@ namespace { } /// Compute the encoding for a set of VolIDs within a readout descriptor - static pair<VolumeID, VolumeID> encoding(const IDDescriptor iddesc, const VolIDs& ids) { + static Encoding update_encoding(const IDDescriptor iddesc, const VolIDs& ids, const Encoding& initial) { + VolumeID volume_id = initial.first, mask = initial.second; + for (VolIDs::const_iterator i = ids.begin(); i != ids.end(); ++i) { + const PlacedVolume::VolID& id = (*i); + IDDescriptor::Field f = iddesc.field(id.first); + VolumeID msk = f->mask(); + int offset = f->offset(); + volume_id |= f->value(id.second << offset) << offset; + mask |= msk; + } + return make_pair(volume_id, mask); + } + /// Compute the encoding for a set of VolIDs within a readout descriptor + static Encoding encoding(const IDDescriptor iddesc, const VolIDs& ids) { VolumeID volume_id = 0, mask = 0; for (VolIDs::const_iterator i = ids.begin(); i != ids.end(); ++i) { const PlacedVolume::VolID& id = (*i); @@ -165,13 +277,47 @@ namespace { } return make_pair(volume_id, mask); } + + void add_entry(SensitiveDetector sd, DetElement parent, DetElement e, + const TGeoNode* n, const Encoding& code, Chain& nodes) + { + if ( sd.isValid() ) { + if (m_entries.find(code.first) == m_entries.end()) { + Readout ro = sd.readout(); + string sd_name = sd.name(); + DetElement sub_detector = m_lcdd.detector(sd_name); + VolumeManager section = m_volManager.addSubdetector(sub_detector, ro); + // This is the block, we effectively have to save for each physical volume with a VolID + VolumeManager::Context* context = new VolumeManager::Context; + context->identifier = code.first; + context->mask = code.second; + context->detector = parent; + context->placement = PlacedVolume(n); + context->element = e; + //context->volID = ids; + context->path = nodes; + for (size_t i = nodes.size(); i > 1; --i) { // Omit the placement of the parent DetElement + TGeoMatrix* m = nodes[i - 1]->GetMatrix(); + context->toWorld.MultiplyLeft(m); + } + context->toDetector = context->toWorld; + context->toDetector.MultiplyLeft(nodes[0]->GetMatrix()); + context->toWorld.MultiplyLeft(&parent.worldTransformation()); + if (!section.adoptPlacement(context)) { + print_node(sd, parent, e, n, code, nodes); + } + m_entries.insert(code.first); + } + } + } + void add_entry(SensitiveDetector sd, DetElement parent, DetElement e, const TGeoNode* n, const VolIDs& ids, Chain& nodes) { if ( sd.isValid() ) { Readout ro = sd.readout(); IDDescriptor iddesc = ro.idSpec(); - pair<VolumeID, VolumeID> code = encoding(iddesc, ids); + Encoding code = encoding(iddesc, ids); if (m_entries.find(code.first) == m_entries.end()) { string sd_name = sd.name(); @@ -209,7 +355,7 @@ namespace { const IDDescriptor& en = ro.idSpec(); PlacedVolume pv = Ref_t(n); bool sensitive = pv.volume().isSensitive(); - pair<VolumeID, VolumeID> code = encoding(en, ids); + Encoding code = encoding(en, ids); VolumeID volume_id = code.first; //if ( !sensitive ) return; @@ -233,6 +379,30 @@ namespace { for(vector<const TGeoNode*>::const_iterator j=nodes.begin(); j!=nodes.end();++j) log << (void*)(*j) << " "; printout(DEBUG,"VolumeManager",log.str().c_str()); +#endif + } + void print_node(SensitiveDetector sd, DetElement parent, DetElement e, + const TGeoNode* n, const Encoding& code, const Chain& /* nodes */) const + { + static int s_count = 0; + Readout ro = sd.readout(); + PlacedVolume pv = Ref_t(n); + bool sensitive = pv.volume().isSensitive(); + VolumeID volume_id = code.first; + + //if ( !sensitive ) return; + ++s_count; + stringstream log; + log << s_count << ": " << parent.name() << ": " << e.name() + << " ro:" << ro.ptr() << " pv:" << n->GetName() << " id:" + << (void*) volume_id << " Sensitive:" << yes_no(sensitive); + printout(DEBUG, "VolumeManager", log.str().c_str()); +#if 0 + log.str(""); + log << s_count << ": " << e.name() << " Detector GeoNodes:"; + for(vector<const TGeoNode*>::const_iterator j=nodes.begin(); j!=nodes.end();++j) + log << (void*)(*j) << " "; + printout(DEBUG,"VolumeManager",log.str().c_str()); #endif } }; diff --git a/DDCore/src/XML/XMLDimension.cpp b/DDCore/src/XML/XMLDimension.cpp index e23ec17d7..b3cad12ce 100644 --- a/DDCore/src/XML/XMLDimension.cpp +++ b/DDCore/src/XML/XMLDimension.cpp @@ -124,6 +124,11 @@ XML_ATTR_ACCESSOR(double, phi_tilt) XML_ATTR_ACCESSOR(int, nphi) XML_ATTR_ACCESSOR(double, rc) XML_ATTR_ACCESSOR(int, nz) +XML_ATTR_ACCESSOR(int, key) +XML_ATTR_ACCESSOR(int, key_min) +XML_ATTR_ACCESSOR(int, key_max) +XML_ATTR_ACCESSOR(int, key_val) +XML_ATTR_ACCESSOR(int, key_value) XML_ATTR_ACCESSOR(double, zstart) XML_ATTR_ACCESSOR_DOUBLE(start) XML_ATTR_ACCESSOR_DOUBLE(end) diff --git a/DDCore/src/plugins/Compact2Objects.cpp b/DDCore/src/plugins/Compact2Objects.cpp index df98f35e0..e308ef7ae 100644 --- a/DDCore/src/plugins/Compact2Objects.cpp +++ b/DDCore/src/plugins/Compact2Objects.cpp @@ -20,6 +20,7 @@ #include "DD4hep/Printout.h" #include "DD4hep/Plugins.h" #include "DD4hep/objects/DetectorInterna.h" +#include "DD4hep/objects/ObjectsInterna.h" #include "XML/DocumentHandler.h" #include "XML/Utilities.h" @@ -57,6 +58,7 @@ namespace DD4hep { template <> void Converter<AlignmentEntry>::operator()(xml_h e) const; template <> void Converter<Region>::operator()(xml_h e) const; template <> void Converter<Readout>::operator()(xml_h e) const; + template <> void Converter<Segmentation>::operator()(xml_h e) const; template <> void Converter<LimitSet>::operator()(xml_h e) const; template <> void Converter<Property>::operator()(xml_h e) const; template <> void Converter<CartesianField>::operator()(xml_h e) const; @@ -533,6 +535,87 @@ template <> void Converter<Region>::operator()(xml_h elt) const { lcdd.addRegion(region); } + +/** Specialized converter for compact readout objects. + * + * <readout name="HcalBarrelHits"> + * <segmentation type="RegularNgonCartesianGridXY" gridSizeX="3.0*cm" gridSizeY="3.0*cm" /> + * <id>system:6,barrel:3,module:4,layer:8,slice:5,x:32:-16,y:-16</id> + * </readout> + */ +template <> void Converter<Segmentation>::operator()(xml_h seg) const { + string type = seg.attr<string>(_U(type)); + string name = seg.hasAttr(_U(name)) ? seg.attr<string>(_U(name)) : string(); + SegmentationObject** object = _option<SegmentationObject*>(); + Segmentation segment(type, name); + if (segment.isValid()) { + typedef Segmentation::Parameters _PARS; + const _PARS& pars = segment.parameters(); + for(_PARS::const_iterator it = pars.begin(); it != pars.end(); ++it) { + Segmentation::Parameter p = *it; + XML::Strng_t pNam(p->name()); + if ( seg.hasAttr(pNam) ) { + string pType = p->type(); + if ( pType.compare("int") == 0 ) { + typedef DDSegmentation::TypedSegmentationParameter<int> ParInt; + static_cast<ParInt*>(p)->setTypedValue(seg.attr<int>(pNam)); + } else if ( pType.compare("float") == 0 ) { + typedef DDSegmentation::TypedSegmentationParameter<float> ParFloat; + static_cast<ParFloat*>(p)->setTypedValue(seg.attr<float>(pNam)); + } else if ( pType.compare("doublevec") == 0 ) { + vector<double> valueVector; + string par = seg.attr<string>(pNam); + printout(DEBUG, "Compact", "++ Converting this string structure: %s.",par.c_str()); + vector<string> elts = DDSegmentation::splitString(par); + for (vector<string>::const_iterator j = elts.begin(); j != elts.end(); ++j) { + if ((*j).empty()) continue; + valueVector.push_back(_toDouble((*j))); + } + typedef DDSegmentation::TypedSegmentationParameter< vector<double> > ParDouVec; + static_cast<ParDouVec*>(p)->setTypedValue(valueVector); + } else if ( pType.compare("double" ) == 0) { + typedef DDSegmentation::TypedSegmentationParameter<double>ParDouble; + static_cast<ParDouble*>(p)->setTypedValue(seg.attr<double>(pNam)); + } else { + p->setValue(seg.attr<string>(pNam)); + } + } else if (not p->isOptional()) { + throw_print("FAILED to create segmentation: " + type + ". Missing mandatory parameter: " + p->name() + "!"); + } + } + long key_min = 0, key_max = 0; + DDSegmentation::Segmentation* base = segment->segmentation; + for(xml_coll_t sub(seg,_U(segmentation)); sub; ++sub) { + SegmentationObject* sub_object = 0; + Converter<Segmentation> sub_conv(lcdd,param,&sub_object); + sub_conv(sub); + if ( sub_object ) { + xml_dim_t s(sub); + if ( sub.hasAttr(_U(key_value)) ) { + key_min = key_max = s.key_value(); + } + else if ( sub.hasAttr(_U(key_min)) && sub.hasAttr(_U(key_max)) ) { + key_min = s.key_min(); + key_max = s.key_max(); + } + else { + stringstream tree; + XML::dump_tree(sub,tree); + throw_print("Nested segmentations: Invalid key specification:"+tree.str()); + } + printout(DEBUG,"Compact","++ Segmentation [%s/%s]: Add sub-segmentation %s [%s]", + name.c_str(), type.c_str(), + sub_object->segmentation->name().c_str(), + sub_object->segmentation->type().c_str()); + base->addSubsegmentation(key_min, key_max, sub_object->segmentation); + sub_object->segmentation = 0; + delete sub_object; + } + } + } + if ( object ) *object = segment.ptr(); +} + /** Specialized converter for compact readout objects. * * <readout name="HcalBarrelHits"> @@ -547,47 +630,14 @@ template <> void Converter<Readout>::operator()(xml_h e) const { printout(DEBUG, "Compact", "++ Converting readout structure: %s.",ro.name()); if (seg) { // Segmentation is not mandatory! - string type = seg.attr<string>(_U(type)); - Segmentation segment(type, name); - if (segment.isValid()) { - typedef Segmentation::Parameters _PARS; - segment->parameters(); - const _PARS& parameters = segment.parameters(); - _PARS::const_iterator it; - for (it = parameters.begin(); it != parameters.end(); ++it) { - Segmentation::Parameter p = *it; - XML::Strng_t pNam(p->name()); - if ( seg.hasAttr(pNam) ) { - string pType = p->type(); - if ( pType.compare("int") == 0 ) { - typedef DDSegmentation::TypedSegmentationParameter<int> ParInt; - static_cast<ParInt*>(p)->setTypedValue(seg.attr<int>(pNam)); - } else if ( pType.compare("float") == 0 ) { - typedef DDSegmentation::TypedSegmentationParameter<float> ParFloat; - static_cast<ParFloat*>(p)->setTypedValue(seg.attr<float>(pNam)); - } else if ( pType.compare("doublevec") == 0 ) { - vector<double> valueVector; - string param = seg.attr<string>(pNam); - printout(DEBUG, "Compact", "++ Converting this string structure: %s.",param.c_str()); - vector<string> elts = DDSegmentation::splitString(param); - for (vector<string>::const_iterator j = elts.begin(); j != elts.end(); ++j) { - if ((*j).empty()) continue; - valueVector.push_back(_toDouble((*j))); - } - typedef DDSegmentation::TypedSegmentationParameter< vector<double> > ParDouVec; - static_cast<ParDouVec*>(p)->setTypedValue(valueVector); - } else if ( pType.compare("double" ) == 0) { - typedef DDSegmentation::TypedSegmentationParameter<double>ParDouble; - static_cast<ParDouble*>(p)->setTypedValue(seg.attr<double>(pNam)); - } else { - p->setValue(seg.attr<string>(pNam)); - } - } else if (not p->isOptional()) { - throw_print("FAILED to create segmentation: " + type + ". Missing mandatory parameter: " + p->name() + "!"); - } - } + SegmentationObject* object = 0; + Converter<Segmentation> converter(lcdd,param,&object); + converter(seg); + if ( object ) { + object->setName(name); + Segmentation segment(object); + ro.setSegmentation(segment); } - ro.setSegmentation(segment); } xml_h id = e.child(_U(id)); if (id) { @@ -597,6 +647,36 @@ template <> void Converter<Readout>::operator()(xml_h e) const { ro.setIDDescriptor(idSpec); lcdd.addIDSpecification(idSpec); } + for(xml_coll_t colls(e,_U(hits_collections)); colls; ++colls) { + string hits_key; + if ( colls.hasAttr(_U(key)) ) hits_key = colls.attr<string>(_U(key)); + for(xml_coll_t coll(colls,_U(hits_collection)); coll; ++coll) { + xml_dim_t c(coll); + string coll_name = c.nameStr(); + string coll_key = hits_key; + long key_min = 0, key_max = 0; + + if ( c.hasAttr(_U(key)) ) { + coll_key = c.attr<string>(_U(key)); + } + if ( c.hasAttr(_U(key_value)) ) { + key_max = key_min = c.key_value(); + } + else if ( c.hasAttr(_U(key_min)) && c.hasAttr(_U(key_max)) ) { + key_min = c.key_min(); + key_max = c.key_max(); + } + else { + stringstream tree; + XML::dump_tree(e,tree); + throw_print("Reaout: Invalid specificatrion for multiple hit collections."+tree.str()); + } + printout(DEBUG,"Compact","++ Readout[%s]: Add hit collection %s [%s] %d-%d", + ro.name(), coll_name.c_str(), coll_key.c_str(), key_min, key_max); + HitCollection hits(coll_name, coll_key, key_min, key_max); + ro->hits.push_back(hits); + } + } lcdd.addReadout(ro); } @@ -823,6 +903,9 @@ template <> void Converter<DetElement>::operator()(xml_h element) const { DetElement det(Ref_t(PluginService::Create<NamedObject*>(type, &lcdd, &element, &sens))); if (det.isValid()) { setChildTitles(make_pair(name, det)); + if ( sd.isValid() ) { + det->flag |= DetElement::Object::HAVE_SENSITIVE_DETECTOR; + } } printout(det.isValid() ? INFO : ERROR, "Compact", "%s subdetector:%s of type %s %s", (det.isValid() ? "++ Converted" : "FAILED "), name.c_str(), type.c_str(), diff --git a/DDDetectors/src/CylindricalEndcapCalorimeter_geo.cpp b/DDDetectors/src/CylindricalEndcapCalorimeter_geo.cpp index 8ceb3201a..b05f5795d 100644 --- a/DDDetectors/src/CylindricalEndcapCalorimeter_geo.cpp +++ b/DDDetectors/src/CylindricalEndcapCalorimeter_geo.cpp @@ -16,6 +16,7 @@ // //========================================================================== #include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/Printout.h" #include "XML/Layering.h" using namespace std; @@ -70,6 +71,7 @@ static Ref_t create_detector(LCDD& lcdd, xml_h e, SensitiveDetector sens) { Position layer_pos(0,0,zlayer-zmin-totWidth/2+layerWidth/2); pv = envelopeVol.placeVolume(layer_vol,layer_pos); pv.addPhysVolID("layer",layer_num); + printout(DEBUG,"Calo","CylindricalEndcapCalorimeter: built layer %d -> %s",layer_num,layer_name.c_str()); ++layer_num; } } diff --git a/DDG4/include/DDG4/Geant4Action.h b/DDG4/include/DDG4/Geant4Action.h index a21b00df9..76dddbc02 100644 --- a/DDG4/include/DDG4/Geant4Action.h +++ b/DDG4/include/DDG4/Geant4Action.h @@ -155,6 +155,9 @@ namespace DD4hep { void add(T* obj) { m_v.push_back(obj); } + void add_front(T* obj) { + m_v.insert(m_v.begin(), obj); + } operator const _V&() const { return m_v; } diff --git a/DDG4/include/DDG4/Geant4ReadoutVolumeFilter.h b/DDG4/include/DDG4/Geant4ReadoutVolumeFilter.h new file mode 100644 index 000000000..bf93d2078 --- /dev/null +++ b/DDG4/include/DDG4/Geant4ReadoutVolumeFilter.h @@ -0,0 +1,61 @@ +// $Id$ +//========================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== +#ifndef DD4HEP_DDG4_GEANT4READOUTVOLUMEFILTER_H +#define DD4HEP_DDG4_GEANT4READOUTVOLUMEFILTER_H + +// Framework include files +#include "DD4hep/Readout.h" +#include "DD4hep/IDDescriptor.h" +#include "DDG4/Geant4SensDetAction.h" + +/// Namespace for the AIDA detector description toolkit +namespace DD4hep { + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace Simulation { + + /// Default base class for all geant 4 tracking actions. + /** + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_SIMULATION + */ + class Geant4ReadoutVolumeFilter: public Geant4Filter { + typedef Geometry::IDDescriptor::Field Field; + typedef Geometry::Readout Readout; + typedef Readout::Collection Collection; + + protected: + /// Reference to readout descriptor + Readout m_readout; + /// Collection index + const Collection* m_collection; + /// Bit field value from ID descriptor + Field m_key; + + public: + /// Standard constructor + Geant4ReadoutVolumeFilter(Geant4Context* context, + const std::string& name, + Readout ro, + const std::string& coll); + /// Default destructor + virtual ~Geant4ReadoutVolumeFilter(); + /// Filter action. Return true if hits should be processed + virtual bool operator()(const G4Step* step) const; + }; + } // End namespace Simulation +} // End namespace DD4hep + +#endif // DD4HEP_DDG4_GEANT4READOUTVOLUMEFILTER_H diff --git a/DDG4/include/DDG4/Geant4SensDetAction.h b/DDG4/include/DDG4/Geant4SensDetAction.h index c1034b29a..60b333071 100644 --- a/DDG4/include/DDG4/Geant4SensDetAction.h +++ b/DDG4/include/DDG4/Geant4SensDetAction.h @@ -207,9 +207,15 @@ namespace DD4hep { /// Add an actor responding to all callbacks. Sequence takes ownership. void adopt(Geant4Filter* filter); + /// Add an actor responding to all callbacks to the sequence front. Sequence takes ownership. + void adopt_front(Geant4Filter* filter); + /// Add an actor responding to all callbacks. Sequence takes ownership. void adoptFilter(Geant4Action* filter); + /// Add an actor responding to all callbacks to the sequence front. Sequence takes ownership. + void adoptFilter_front(Geant4Action* filter); + /// Callback before hit processing starts. Invoke all filters. /** Return fals if any filter returns false */ @@ -227,6 +233,9 @@ namespace DD4hep { /// Retrieve the hits collection associated with this detector by its collection identifier HitCollection* collectionByID(size_t id); + /// Define collections created by this sensitivie action object + virtual void defineCollections(); + /// G4VSensitiveDetector interface: Method invoked at the begining of each event. /** The hits collection(s) created by this sensitive detector must * be set to the G4HCofThisEvent object at one of these two methods. @@ -460,10 +469,12 @@ namespace DD4hep { public: typedef T UserData; protected: + /// Property: collection name. If not set default is readout name! + std::string m_collectionName; /// Collection identifier - size_t m_collectionID; + size_t m_collectionID; /// User data block - UserData m_userData; + UserData m_userData; public: /// Standard , initializing constructor @@ -473,12 +484,31 @@ namespace DD4hep { Geometry::LCDD& lcdd); /// Default destructor virtual ~Geant4SensitiveAction(); + + /// Define collections created by this sensitivie action object + virtual void defineCollections(); + + /// Define readout specific hit collection with volume ID filtering + /** + * Convenience function. To be called by specialized sensitive actions inheriting this class. + * + * - If the property CollectionName ist NOT set, one default collection is declared, + * with the readout name as the collection name. + * - If the property CollectionName is set, the readout structure is searched for + * the corresponding entry and a collection with the corrsponding name is declared. + * At the same time a VolumeID filter is injected at the front of the sensitive's + * filter queue to ONLY act on volume IDs matching this criterium. + */ + template <typename HIT> size_t declareReadoutFilteredCollection(); + + /// Define readout specific hit collection. matching name must be present in readout structure + template <typename HIT> + size_t defineReadoutCollection(const std::string collection_name); + /// Initialization overload for specialization virtual void initialize(); /// Finalization overload for specialization virtual void finalize(); - /// Define collections created by this sensitivie action object - virtual void defineCollections() {} /// G4VSensitiveDetector interface: Method invoked at the begining of each event. virtual void begin(G4HCofThisEvent* hce); /// G4VSensitiveDetector interface: Method invoked at the end of each event. diff --git a/DDG4/include/DDG4/Geant4SensDetAction.inl b/DDG4/include/DDG4/Geant4SensDetAction.inl index bdf3e6aa0..cb78c8ac9 100644 --- a/DDG4/include/DDG4/Geant4SensDetAction.inl +++ b/DDG4/include/DDG4/Geant4SensDetAction.inl @@ -16,10 +16,16 @@ #ifndef DD4HEP_DDG4_GEANT4SENSDETACTION_H #include "DDG4/Geant4SensDetAction.h" #endif +#include "DD4hep/objects/ObjectsInterna.h" #include "DD4hep/InstanceCount.h" + +#include "DDG4/Geant4ReadoutVolumeFilter.h" #include "DDG4/Geant4Data.h" +/// Namespace for the AIDA detector description toolkit namespace DD4hep { + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit namespace Simulation { /// Standard , initializing constructor @@ -28,10 +34,10 @@ namespace DD4hep { const std::string& nam, Geometry::DetElement det, Geometry::LCDD& lcdd_ref) - : Geant4Sensitive(ctxt,nam,det,lcdd_ref), m_collectionID(0) + : Geant4Sensitive(ctxt,nam,det,lcdd_ref), m_collectionName(), m_collectionID(0) { + declareProperty("CollectionName", m_collectionName); initialize(); - defineCollections(); InstanceCount::increment(this); } @@ -69,10 +75,47 @@ namespace DD4hep { Geant4Sensitive::clear(hce); } + /// Define collections created by this sensitivie action object + template <typename T> void Geant4SensitiveAction<T>::defineCollections() { + } + + /// Define readout specific hit collection. matching name must be present in readout structure + template <typename T> template <typename HIT> + size_t Geant4SensitiveAction<T>::defineReadoutCollection(const std::string match_name) { + typedef Geometry::ReadoutObject _RO; + Readout ro = m_sensitive.readout(); + for(_RO::Collections::const_iterator i=ro->hits.begin(); i!=ro->hits.end(); ++i) { + const _RO::Collection& c = *i; + if ( c.name == match_name ) { + size_t coll_id = defineCollection<HIT>(match_name); + Geant4Filter* filter = + new Geant4ReadoutVolumeFilter(context(),name()+"_"+c.name,ro,c.name); + adopt_front(filter); + return coll_id; + } + } + except("+++ Custom collection name '%s' not defined in the Readout object: %s.", + m_collectionName.c_str(), ro.name()); + return 0; // Anyhow not reachable. Exception thrown before + } + + /// Define readout specific hit collection with volume ID filtering + template <typename T> template <typename HIT> + size_t Geant4SensitiveAction<T>::declareReadoutFilteredCollection() + { + if ( m_collectionName.empty() ) { + Readout ro = m_sensitive.readout(); + return defineCollection<HIT>(ro.name()); + } + return defineReadoutCollection<HIT>(m_collectionName); + } + // Forward declarations typedef Geant4HitData::Contribution HitContribution; - } -} + + + } // End namespace Simulation +} // End namespace DD4hep #include "DD4hep/Printout.h" diff --git a/DDG4/plugins/Geant4SDActions.cpp b/DDG4/plugins/Geant4SDActions.cpp index 1f7532115..f06639fa9 100644 --- a/DDG4/plugins/Geant4SDActions.cpp +++ b/DDG4/plugins/Geant4SDActions.cpp @@ -41,8 +41,8 @@ namespace DD4hep { */ /// Define collections created by this sensitivie action object - template <> void Geant4SensitiveAction<Geant4Tracker>::defineCollections() { - m_collectionID = defineCollection<Geant4Tracker::Hit>(m_sensitive.readout().name()); + template <> void Geant4SensitiveAction<Geant4Tracker>::defineCollections() { + m_collectionID = declareReadoutFilteredCollection<Geant4Tracker::Hit>(); } /// Method for generating hit(s) using the information of G4Step object. @@ -98,8 +98,9 @@ namespace DD4hep { /// Define collections created by this sensitivie action object template <> void Geant4SensitiveAction<Geant4Calorimeter>::defineCollections() { - m_collectionID = defineCollection<Geant4Calorimeter::Hit>(m_sensitive.readout().name()); + m_collectionID = declareReadoutFilteredCollection<Geant4Calorimeter::Hit>(); } + /// Method for generating hit(s) using the information of G4Step object. template <> bool Geant4SensitiveAction<Geant4Calorimeter>::process(G4Step* step,G4TouchableHistory*) { typedef Geant4Calorimeter::Hit Hit; @@ -135,8 +136,9 @@ namespace DD4hep { hit = new Hit(global); hit->cellID = cell; coll->add(hit); - printM2("%s> CREATE hit with deposit:%e MeV Pos:%8.2f %8.2f %8.2f %s", - c_name(),contrib.deposit,pos.X,pos.Y,pos.Z,handler.path().c_str()); + printM2("%s> CREATE hit with deposit:%e MeV Pos:%8.2f %8.2f %8.2f %s [%s]", + c_name(),contrib.deposit,pos.X,pos.Y,pos.Z,handler.path().c_str(), + coll->GetName().c_str()); if ( 0 == hit->cellID ) { // for debugging only! hit->cellID = cellID(step); except("+++ Invalid CELL ID for hit!"); @@ -173,7 +175,7 @@ namespace DD4hep { /// Define collections created by this sensitivie action object template <> void Geant4SensitiveAction<Geant4OpticalCalorimeter>::defineCollections() { - m_collectionID = defineCollection<Geant4Calorimeter::Hit>(m_sensitive.readout().name()); + m_collectionID = declareReadoutFilteredCollection<Geant4Calorimeter::Hit>(); } /// Method for generating hit(s) using the information of G4Step object. template <> bool Geant4SensitiveAction<Geant4OpticalCalorimeter>::process(G4Step* step,G4TouchableHistory*) { @@ -235,7 +237,7 @@ namespace DD4hep { /// Define collections created by this sensitivie action object template <> void Geant4SensitiveAction<Geant4ScintillatorCalorimeter>::defineCollections() { - m_collectionID = defineCollection<Geant4Calorimeter::Hit>(m_sensitive.readout().name()); + m_collectionID = declareReadoutFilteredCollection<Geant4Calorimeter::Hit>(); } /// Method for generating hit(s) using the information of G4Step object. template <> bool Geant4SensitiveAction<Geant4ScintillatorCalorimeter>::process(G4Step* step,G4TouchableHistory*) { @@ -460,7 +462,7 @@ namespace DD4hep { /// Define collections created by this sensitivie action object template <> void Geant4SensitiveAction<TrackerCombine>::defineCollections() { - m_collectionID = defineCollection<Geant4Tracker::Hit>(m_sensitive.readout().name()); + m_collectionID = declareReadoutFilteredCollection<Geant4Tracker::Hit>(); } /// Method for generating hit(s) using the information of G4Step object. diff --git a/DDG4/plugins/Geant4TrackerWeightedSD.cpp b/DDG4/plugins/Geant4TrackerWeightedSD.cpp index 787736116..80ff2d75a 100644 --- a/DDG4/plugins/Geant4TrackerWeightedSD.cpp +++ b/DDG4/plugins/Geant4TrackerWeightedSD.cpp @@ -396,7 +396,7 @@ namespace DD4hep { /// Define collections created by this sensitivie action object template <> void Geant4SensitiveAction<TrackerWeighted>::defineCollections() { - m_collectionID = defineCollection<Geant4Tracker::Hit>(m_sensitive.readout().name()); + m_collectionID = declareReadoutFilteredCollection<Geant4Tracker::Hit>(); } /// Method for generating hit(s) using the information of G4Step object. diff --git a/DDG4/python/DDG4.py b/DDG4/python/DDG4.py index c34b82d7a..61ce223a4 100644 --- a/DDG4/python/DDG4.py +++ b/DDG4/python/DDG4.py @@ -484,36 +484,79 @@ class Geant4: sdtyp = self.sensitive_types[typ] print '+++ %-32s type:%-12s --> Sensitive type: %s'%(o.name(), typ, sdtyp,) - def setupDetector(self,name,action): + def setupSensitiveSequencer(self, name, action): + if isinstance( action, tuple ): + sensitive_type = action[0] + else: + sensitive_type = action + seq = SensitiveSequence(self.kernel(),'Geant4SensDetActionSequence/'+name) + seq.enableUI() + return seq + + def setupDetector(self,name,action,collections=None): #fg: allow the action to be a tuple with parameter dictionary sensitive_type = "" parameterDict = {} - if isinstance( action, tuple ): + if isinstance(action,tuple) or isinstance(action,list): sensitive_type = action[0] parameterDict = action[1] else: sensitive_type = action - + seq = SensitiveSequence(self.kernel(),'Geant4SensDetActionSequence/'+name) - act = SensitiveAction(self.kernel(),sensitive_type+'/'+name+'Handler',name) seq.enableUI() - act.enableUI() - seq.add(act) - for parameter, value in parameterDict.iteritems(): - setattr( act, parameter, value) - return (seq,act) - - def setupCalorimeter(self,name,type=None): + acts = [] + if collections is None: + sd = self.lcdd.sensitiveDetector(name) + ro = sd.readout() + #print dir(ro) + collections = ro.collectionNames() + if len(collections)==0: + act = SensitiveAction(self.kernel(),sensitive_type+'/'+name+'Handler',name) + for parameter, value in parameterDict.iteritems(): + setattr( act, parameter, value) + acts.append(act) + + # Work down the collections if present + if collections is not None: + for coll in collections: + params = {} + if isinstance(coll,tuple) or isinstance(coll,list): + if len(coll)>2: + coll_nam = coll[0] + sensitive_type = coll[1] + params = coll[2] + elif len(coll)>1: + coll_nam = coll[0] + sensitive_type = coll[1] + else: + coll_nam = coll[0] + else: + coll_nam = coll + act = SensitiveAction(self.kernel(),sensitive_type+'/'+coll_nam+'Handler',name) + act.CollectionName = coll_nam + for parameter, value in params.iteritems(): + setattr( act, parameter, value) + acts.append(act) + + for act in acts: + act.enableUI() + seq.add(act) + if len(acts)>1: + return (seq,acts) + return (seq,acts[0]) + + def setupCalorimeter(self,name,type=None,collections=None): sd = self.lcdd.sensitiveDetector(name) sd.setType('calorimeter') if type is None: type = self.sensitive_types['calorimeter'] - return self.setupDetector(name,type) + return self.setupDetector(name,type,collections) - def setupTracker(self,name,type=None): + def setupTracker(self,name,type=None,collections=None): sd = self.lcdd.sensitiveDetector(name) sd.setType('tracker') if type is None: type = self.sensitive_types['tracker'] - return self.setupDetector(name,type) + return self.setupDetector(name,type,collections) def setupTrackingField(self, name='MagFieldTrackingSetup', stepper='HelixSimpleRunge', equation='Mag_UsualEqRhs',prt=False): import SystemOfUnits diff --git a/DDG4/src/Geant4OutputAction.cpp b/DDG4/src/Geant4OutputAction.cpp index 016afbc1e..ea64cfdf2 100644 --- a/DDG4/src/Geant4OutputAction.cpp +++ b/DDG4/src/Geant4OutputAction.cpp @@ -34,8 +34,8 @@ Geant4OutputAction::Geant4OutputAction(Geant4Context* ctxt, const string& nam) InstanceCount::increment(this); declareProperty("Output", m_output); declareProperty("HandleErrorsAsFatal", m_errorFatal=true); - //ctxt->runAction().callAtBegin(this, &Geant4OutputAction::beginRun); - //ctxt->runAction().callAtEnd(this, &Geant4OutputAction::endRun); + // Need to instantiate run action to configure fibers + if ( 0 == &ctxt->runAction() ) {} } /// Default destructor diff --git a/DDG4/src/Geant4ReadoutVolumeFilter.cpp b/DDG4/src/Geant4ReadoutVolumeFilter.cpp new file mode 100644 index 000000000..f6aff8b1b --- /dev/null +++ b/DDG4/src/Geant4ReadoutVolumeFilter.cpp @@ -0,0 +1,60 @@ +// $Id$ +//========================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== + +// Framework include files +#include "DD4hep/Readout.h" +#include "DD4hep/InstanceCount.h" +#include "DD4hep/objects/ObjectsInterna.h" +#include "DDG4/Geant4Mapping.h" +#include "DDG4/Geant4StepHandler.h" +#include "DDG4/Geant4VolumeManager.h" +#include "DDG4/Geant4ReadoutVolumeFilter.h" + +using namespace DD4hep::Simulation; + +/// Standard constructor +Geant4ReadoutVolumeFilter::Geant4ReadoutVolumeFilter(Geant4Context* ctxt, + const std::string& nam, + Readout ro, + const std::string& coll) + : Geant4Filter(ctxt, nam), m_readout(ro), m_collection(0), m_key(0) +{ + InstanceCount::increment(this); + for(size_t i=0; i<ro->hits.size(); ++i) { + const Collection& c = ro->hits[i]; + if ( c.name == coll ) { + m_collection = &c; + m_key = ro.idSpec().field(c.key); + return; + } + } + except("+++ Custom collection name '%s' not defined in the Readout object: %s.", + coll.c_str(), ro.name()); +} + +/// Default destructor +Geant4ReadoutVolumeFilter::~Geant4ReadoutVolumeFilter() { + InstanceCount::decrement(this); +} + +/// Filter action. Return true if hits should be processed +bool Geant4ReadoutVolumeFilter::operator()(const G4Step* s) const { + Geant4StepHandler step(s); + Geant4VolumeManager volMgr = Geant4Mapping::instance().volumeManager(); + VolumeID id = volMgr.volumeID(step.preTouchable()); + long64 key = m_key->value(id); + if ( m_collection->key_min <= key && m_collection->key_max >= key ) + return true; + return false; +} diff --git a/DDG4/src/Geant4SensDetAction.cpp b/DDG4/src/Geant4SensDetAction.cpp index 05552e435..8d431d65c 100644 --- a/DDG4/src/Geant4SensDetAction.cpp +++ b/DDG4/src/Geant4SensDetAction.cpp @@ -121,7 +121,23 @@ void Geant4Sensitive::adopt(Geant4Filter* filter) { m_filters.add(filter); return; } - throw runtime_error("Geant4SensDetActionSequence: Attempt to add invalid sensitive filter!"); + throw runtime_error("Geant4Sensitive: Attempt to add invalid sensitive filter!"); +} + +/// Add an actor responding to all callbacks. Sequence takes ownership. +void Geant4Sensitive::adoptFilter_front(Geant4Action* action) { + Geant4Filter* filter = dynamic_cast<Geant4Filter*>(action); + adopt_front(filter); +} + +/// Add an actor responding to all callbacks. Sequence takes ownership. +void Geant4Sensitive::adopt_front(Geant4Filter* filter) { + if (filter) { + filter->addRef(); + m_filters.add_front(filter); + return; + } + throw runtime_error("Geant4Sensitive: Attempt to add invalid sensitive filter!"); } /// Callback before hit processing starts. Invoke all filters. @@ -165,6 +181,10 @@ Geant4HitCollection* Geant4Sensitive::collectionByID(size_t id) { return sequence().collectionByID(id); } +/// Define collections created by this sensitivie action object +void Geant4Sensitive::defineCollections() { +} + /// Method invoked at the begining of each event. void Geant4Sensitive::begin(G4HCofThisEvent* /* HCE */) { } @@ -284,6 +304,7 @@ size_t Geant4SensDetActionSequence::Geant4SensDetActionSequence::defineCollectio size_t count = 0; m_detector = sens_det; m_actors(&Geant4Sensitive::setDetector, sens_det); + m_actors(&Geant4Sensitive::defineCollections); for (HitCollections::const_iterator i = m_collections.begin(); i != m_collections.end(); ++i) { sens_det->defineCollection((*i).first); ++count; diff --git a/DDSegmentation/include/DDSegmentation/MultiSegmentation.h b/DDSegmentation/include/DDSegmentation/MultiSegmentation.h new file mode 100644 index 000000000..557919d89 --- /dev/null +++ b/DDSegmentation/include/DDSegmentation/MultiSegmentation.h @@ -0,0 +1,90 @@ +// $Id$ +//========================================================================== +// AIDA Detector description implementation for LCD +//-------------------------------------------------------------------------- +// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN) +// All rights reserved. +// +// For the licensing terms see $DD4hepINSTALL/LICENSE. +// For the list of contributors see $DD4hepINSTALL/doc/CREDITS. +// +// Author : M.Frank +// +//========================================================================== +#ifndef DDSegmentation_MULITSEGMENTATION_H_ +#define DDSegmentation_MULITSEGMENTATION_H_ + +#include "DDSegmentation/Segmentation.h" + +/// Main handle class to hold a TGeo alignment object of type TGeoPhysicalNode +namespace DD4hep { + + namespace DDSegmentation { + + /// Wrapper to support multiple segmentations + /** + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_DDSEGMENTATION + */ + class MultiSegmentation : public Segmentation { + public: + /// Default constructor passing the encoding string + MultiSegmentation(const std::string& cellEncoding = ""); + + /// Default destructor + virtual ~MultiSegmentation(); + + /// Add subsegmentation. + virtual void addSubsegmentation(long key_min, long key_max, Segmentation* entry); + + /// determine the position based on the cell ID + virtual Vector3D position(const CellID& cellID) const; + + /// determine the cell ID based on the position + virtual CellID cellID(const Vector3D& localPosition, const Vector3D& globalPosition, const VolumeID& volumeID) const; + + /** \brief Returns a vector<double> of the cellDimensions of the given cell ID + in natural order of dimensions, e.g., dx/dy/dz, or dr/r*dPhi + + \param cellID cellID of the cell for which parameters are returned + \return vector<double> in natural order of dimensions, e.g., dx/dy/dz, or dr/r*dPhi + */ + virtual std::vector<double> cellDimensions(const CellID& cellID) const; + + /// access the field name used to discriminate sub-segmentations + const std::string& discriminatorName() const { return m_discriminatorId; } + + /// Discriminating bitfield entry + BitFieldValue* discriminator() const { return m_discriminator; } + + /// Set the underlying decoder + virtual void setDecoder(BitField64* decoder); + + protected: + struct Entry { + long key_min, key_max; + Segmentation* segmentation; + }; + typedef std::vector<Entry> Segmentations; + + /// Access subsegmentation by cell identifier + const Segmentation& subsegmentation(const CellID& cellID) const; + + /// Sub-segmentaion container + Segmentations m_segmentations; + + /// the field name used to discriminate sub-segmentations + std::string m_discriminatorId; + + /// Bitfield corresponding to dicriminator identifier + BitFieldValue* m_discriminator; + + /// Debug flags + int m_debug; + }; + + } /* namespace DDSegmentation */ +} /* namespace DD4hep */ + +#endif /* DDSegmentation_MULITSEGMENTATION_H_ */ diff --git a/DDSegmentation/include/DDSegmentation/Segmentation.h b/DDSegmentation/include/DDSegmentation/Segmentation.h index 68a7bf13d..1db44adca 100644 --- a/DDSegmentation/include/DDSegmentation/Segmentation.h +++ b/DDSegmentation/include/DDSegmentation/Segmentation.h @@ -70,6 +70,8 @@ public: /// Destructor virtual ~Segmentation(); + /// Add subsegmentation. Call only valid for Multi-segmentations. Default implementation throws an exception + virtual void addSubsegmentation(long key_min, long key_max, Segmentation* entry); /// Determine the local position based on the cell ID virtual Vector3D position(const CellID& cellID) const = 0; /// Determine the cell ID based on the position @@ -141,7 +143,7 @@ protected: static int positionToBin(double position, double cellSize, double offset = 0.); /// Helper method to convert a bin number to a 1D position given a vector of binBoundaries - static double binToPosition(CellID bin, std::vector<double> const& cellBoundaries, double offset = 0.); + static double binToPosition(CellID bin, std::vector<double> const& cellBoundaries, double offset = 0.); /// Helper method to convert a 1D position to a cell ID given a vector of binBoundaries static int positionToBin(double position, std::vector<double> const& cellBoundaries, double offset = 0.); diff --git a/DDSegmentation/src/MultiSegmentation.cpp b/DDSegmentation/src/MultiSegmentation.cpp new file mode 100644 index 000000000..1a9d43ec6 --- /dev/null +++ b/DDSegmentation/src/MultiSegmentation.cpp @@ -0,0 +1,91 @@ +/* + * MultiSegmentation.cpp + * + * Created on: Jun 28, 2013 + * Author: Christian Grefe, CERN + */ + +#include "DDSegmentation/MultiSegmentation.h" +#include <iomanip> + +using namespace std; + +namespace DD4hep { + namespace DDSegmentation { + + /// default constructor using an encoding string + MultiSegmentation::MultiSegmentation(const string& cellEncoding) + : Segmentation(cellEncoding), m_discriminator(0), m_debug(0) + { + // define type and description + _type = "MultiSegmentation"; + _description = "Multi-segmenation wrapper segmentation"; + //registerParameter<int>("debug", "Debug flag", m_debug, 0); + registerParameter<string>("key", "Diskriminating field", m_discriminatorId, ""); + } + + /// destructor + MultiSegmentation::~MultiSegmentation() { + for(Segmentations::iterator i=m_segmentations.begin(); i!=m_segmentations.end(); ++i) + delete (*i).segmentation; + m_segmentations.clear(); + } + + /// Add subsegmentation. Call only valid for Multi-segmentations. Default implementation throws an exception + void MultiSegmentation::addSubsegmentation(long key_min, long key_max, Segmentation* entry) { + Entry e; + e.key_min = key_min; + e.key_max = key_max; + e.segmentation = entry; + m_segmentations.push_back(e); + } + + /// Set the underlying decoder + void MultiSegmentation::setDecoder(BitField64* newDecoder) { + this->Segmentation::setDecoder(newDecoder); + for(Segmentations::iterator i=m_segmentations.begin(); i != m_segmentations.end(); ++i) + (*i).segmentation->setDecoder(newDecoder); + m_discriminator = &((*_decoder)[m_discriminatorId]); + } + + /// + const Segmentation& MultiSegmentation::subsegmentation(const CellID& cID) const { + if ( m_discriminator ) { + long seg_id = m_discriminator->value(cID); + for(Segmentations::const_iterator i=m_segmentations.begin(); i != m_segmentations.end(); ++i) { + const Entry& e = *i; + if ( e.key_min<= seg_id && e.key_max >= seg_id ) { + Segmentation* s = e.segmentation; + if ( m_debug > 0 ) { + cout << "MultiSegmentation: id:" << setw(4) << hex << seg_id << dec << " " << s->name(); + const Parameters& pars = s->parameters(); + for(Parameters::const_iterator j=pars.begin(); j!=pars.end();++j) { + cout << " " << (*j)->name() << "=" << (*j)->value(); + } + cout << endl; + } + return *s; + } + } + } + throw runtime_error("MultiSegmentation: Invalid sub-segmentation identifier!");; + } + + /// determine the position based on the cell ID + Vector3D MultiSegmentation::position(const CellID& cID) const { + return subsegmentation(cID).position(cID); + } + + /// determine the cell ID based on the position + CellID MultiSegmentation::cellID(const Vector3D& localPosition, const Vector3D& globalPosition, const VolumeID& vID) const { + return subsegmentation(vID).cellID(localPosition, globalPosition, vID); + } + + vector<double> MultiSegmentation::cellDimensions(const CellID& cID) const { + return subsegmentation(cID).cellDimensions(cID); + } + + REGISTER_SEGMENTATION(MultiSegmentation) + + } /* namespace DDSegmentation */ +} /* namespace DD4hep */ diff --git a/DDSegmentation/src/Segmentation.cpp b/DDSegmentation/src/Segmentation.cpp index 27be7fc59..e9b037f8e 100644 --- a/DDSegmentation/src/Segmentation.cpp +++ b/DDSegmentation/src/Segmentation.cpp @@ -15,185 +15,190 @@ #include <iomanip> namespace DD4hep { -namespace DDSegmentation { - -using std::cerr; -using std::endl; -using std::map; -using std::runtime_error; -using std::stringstream; -using std::vector; - -/// Default constructor used by derived classes passing the encoding string -Segmentation::Segmentation(const std::string& cellEncoding) : - _name("Segmentation"), _type("Segmentation"), _decoder(new BitField64(cellEncoding)), _ownsDecoder(true) { - -} - -/// Default constructor used by derived classes passing an existing decoder -Segmentation::Segmentation(BitField64* newDecoder) : - _name("Segmentation"), _type("Segmentation"), _decoder(newDecoder), _ownsDecoder(false) { -} - -/// Destructor -Segmentation::~Segmentation() { - if (_ownsDecoder and _decoder != 0) { - delete _decoder; - } - map<std::string, SegmentationParameter*>::iterator it; - for (it = _parameters.begin(); it != _parameters.end(); ++it) { - SegmentationParameter* p = it->second; - if (p) { - delete p; - p = 0; - } - } -} - -/// Determine the volume ID from the full cell ID by removing all local fields -VolumeID Segmentation::volumeID(const CellID& cID) const { - map<std::string, StringParameter>::const_iterator it; - _decoder->setValue(cID); - for (it = _indexIdentifiers.begin(); it != _indexIdentifiers.end(); ++it) { - std::string identifier = it->second->typedValue(); - (*_decoder)[identifier] = 0; - } - return _decoder->getValue(); -} - -/// Calculates the neighbours of the given cell ID and adds them to the list of neighbours -void Segmentation::neighbours(const CellID& cID, std::set<CellID>& cellNeighbours) const { - map<std::string, StringParameter>::const_iterator it; - for (it = _indexIdentifiers.begin(); it != _indexIdentifiers.end(); ++it) { - const std::string& identifier = it->second->typedValue(); - _decoder->setValue(cID); - int currentValue = (*_decoder)[identifier]; - // add both neighbouring cell IDs, don't add out of bound indices - try { - (*_decoder)[identifier] = currentValue - 1; - cellNeighbours.insert(_decoder->getValue()); - } catch (runtime_error& e) { - // nothing to do - } - try { - (*_decoder)[identifier] = currentValue + 1; - cellNeighbours.insert(_decoder->getValue()); - } catch (runtime_error& e) { - // nothing to do - } - } -} - -/// Set the underlying decoder -void Segmentation::setDecoder(BitField64* newDecoder) { - if ( _decoder == newDecoder ) return; //self assignment - if (_ownsDecoder) { - delete _decoder; - } - _decoder = newDecoder; - _ownsDecoder = false; -} - -/// Access to parameter by name -Parameter Segmentation::parameter(const std::string& parameterName) const { - map<std::string, Parameter>::const_iterator it; - it = _parameters.find(parameterName); - if (it != _parameters.end()) { - return it->second; - } - stringstream s; - s << "Unknown parameter " << parameterName << " for segmentation type " << _type; - throw std::runtime_error(s.str()); -} - -/// Access to all parameters -Parameters Segmentation::parameters() const { - Parameters pars; - map<std::string, Parameter>::const_iterator it; - for (it = _parameters.begin(); it != _parameters.end(); ++it) { - pars.push_back(it->second); - } - return pars; -} - -/// Set all parameters from an existing set of parameters -void Segmentation::setParameters(const Parameters& pars) { - Parameters::const_iterator it; - for (it = pars.begin(); it != pars.end(); ++it) { - Parameter p = *it; - parameter(p->name())->value() = p->value(); - } -} - -/// Add a cell identifier to this segmentation. Used by derived classes to define their required identifiers -void Segmentation::registerIdentifier(const std::string& idName, const std::string& idDescription, std::string& identifier, - const std::string& defaultValue) { - StringParameter idParameter = - new TypedSegmentationParameter<std::string>(idName, idDescription, identifier, defaultValue, - SegmentationParameter::NoUnit, true); - _parameters[idName] = idParameter; - _indexIdentifiers[idName] = idParameter; -} - -/// Helper method to convert a bin number to a 1D position -double Segmentation::binToPosition(long64 bin, double cellSize, double offset) { - return bin * cellSize + offset; -} - -/// Helper method to convert a 1D position to a cell ID -int Segmentation::positionToBin(double position, double cellSize, double offset) { - if (cellSize <= 1e-10) { - throw runtime_error("Invalid cell size: 0.0"); - } - return int(floor((position + 0.5 * cellSize - offset) / cellSize)); -} - -/// Helper method to convert a bin number to a 1D position given a vector of binBoundaries -double Segmentation::binToPosition(CellID bin, std::vector<double> const& cellBoundaries, double offset) { - return (cellBoundaries[bin+1] + cellBoundaries[bin])*0.5 + offset; -} -/// Helper method to convert a 1D position to a cell ID given a vector of binBoundaries -int Segmentation::positionToBin(double position, std::vector<double> const& cellBoundaries, double offset) { - - // include the lower edge to the segmentation, deal with numerical issues - if(fabs(position/cellBoundaries.front()-1.0) < 3e-12) return 0; - - // include the upper edge of the last bin to the segmentation, deal with numerical issues - if(fabs(position/cellBoundaries.back()-1.0) < 3e-12) return int(cellBoundaries.size()-2); - - // hits outside cannot be treated - if(position < cellBoundaries.front()) { - std::stringstream err; - err << std::setprecision(20) << std::scientific; - err << "Hit Position (" << position << ") is below the acceptance" - << " (" << cellBoundaries.front() << ") " - << "of the segmentation"; - throw std::runtime_error(err.str()); - } - if(position > cellBoundaries.back() ) { - std::stringstream err; - err << std::setprecision(20) << std::scientific; - err << "Hit Position (" << position << ") is above the acceptance" - << " (" << cellBoundaries.back() << ") " - << "of the segmentation"; - throw std::runtime_error(err.str()); - } - - - std::vector<double>::const_iterator bin = std::upper_bound(cellBoundaries.begin(), - cellBoundaries.end(), - position-offset); - - // need to reduce found bin by one, because upper_bound works that way, lower_bound would not help - return bin - cellBoundaries.begin() - 1 ; - -} - -std::vector<double> Segmentation::cellDimensions(const CellID&) const { - std::stringstream errorMessage; - errorMessage << __func__ << " is not implemented for " << _type; - throw std::logic_error(errorMessage.str()); -} - -} /* namespace DDSegmentation */ + namespace DDSegmentation { + + using std::cerr; + using std::endl; + using std::map; + using std::runtime_error; + using std::stringstream; + using std::vector; + + /// Default constructor used by derived classes passing the encoding string + Segmentation::Segmentation(const std::string& cellEncoding) : + _name("Segmentation"), _type("Segmentation"), _decoder(new BitField64(cellEncoding)), _ownsDecoder(true) { + + } + + /// Default constructor used by derived classes passing an existing decoder + Segmentation::Segmentation(BitField64* newDecoder) : + _name("Segmentation"), _type("Segmentation"), _decoder(newDecoder), _ownsDecoder(false) { + } + + /// Destructor + Segmentation::~Segmentation() { + if (_ownsDecoder and _decoder != 0) { + delete _decoder; + } + map<std::string, SegmentationParameter*>::iterator it; + for (it = _parameters.begin(); it != _parameters.end(); ++it) { + SegmentationParameter* p = it->second; + if (p) { + delete p; + p = 0; + } + } + } + + /// Add subsegmentation. Call only valid for Multi-segmentations. Default implementation throws an exception + void Segmentation::addSubsegmentation(long /* key_min */, long /* key_max */, Segmentation* /* entry */) { + throw std::runtime_error("This segmentation type:"+_type+" does not support sub-segmentations."); + } + + /// Determine the volume ID from the full cell ID by removing all local fields + VolumeID Segmentation::volumeID(const CellID& cID) const { + map<std::string, StringParameter>::const_iterator it; + _decoder->setValue(cID); + for (it = _indexIdentifiers.begin(); it != _indexIdentifiers.end(); ++it) { + std::string identifier = it->second->typedValue(); + (*_decoder)[identifier] = 0; + } + return _decoder->getValue(); + } + + /// Calculates the neighbours of the given cell ID and adds them to the list of neighbours + void Segmentation::neighbours(const CellID& cID, std::set<CellID>& cellNeighbours) const { + map<std::string, StringParameter>::const_iterator it; + for (it = _indexIdentifiers.begin(); it != _indexIdentifiers.end(); ++it) { + const std::string& identifier = it->second->typedValue(); + _decoder->setValue(cID); + int currentValue = (*_decoder)[identifier]; + // add both neighbouring cell IDs, don't add out of bound indices + try { + (*_decoder)[identifier] = currentValue - 1; + cellNeighbours.insert(_decoder->getValue()); + } catch (runtime_error& e) { + // nothing to do + } + try { + (*_decoder)[identifier] = currentValue + 1; + cellNeighbours.insert(_decoder->getValue()); + } catch (runtime_error& e) { + // nothing to do + } + } + } + + /// Set the underlying decoder + void Segmentation::setDecoder(BitField64* newDecoder) { + if ( _decoder == newDecoder ) return; //self assignment + if (_ownsDecoder) { + delete _decoder; + } + _decoder = newDecoder; + _ownsDecoder = false; + } + + /// Access to parameter by name + Parameter Segmentation::parameter(const std::string& parameterName) const { + map<std::string, Parameter>::const_iterator it; + it = _parameters.find(parameterName); + if (it != _parameters.end()) { + return it->second; + } + stringstream s; + s << "Unknown parameter " << parameterName << " for segmentation type " << _type; + throw std::runtime_error(s.str()); + } + + /// Access to all parameters + Parameters Segmentation::parameters() const { + Parameters pars; + map<std::string, Parameter>::const_iterator it; + for (it = _parameters.begin(); it != _parameters.end(); ++it) { + pars.push_back(it->second); + } + return pars; + } + + /// Set all parameters from an existing set of parameters + void Segmentation::setParameters(const Parameters& pars) { + Parameters::const_iterator it; + for (it = pars.begin(); it != pars.end(); ++it) { + Parameter p = *it; + parameter(p->name())->value() = p->value(); + } + } + + /// Add a cell identifier to this segmentation. Used by derived classes to define their required identifiers + void Segmentation::registerIdentifier(const std::string& idName, const std::string& idDescription, std::string& identifier, + const std::string& defaultValue) { + StringParameter idParameter = + new TypedSegmentationParameter<std::string>(idName, idDescription, identifier, defaultValue, + SegmentationParameter::NoUnit, true); + _parameters[idName] = idParameter; + _indexIdentifiers[idName] = idParameter; + } + + /// Helper method to convert a bin number to a 1D position + double Segmentation::binToPosition(long64 bin, double cellSize, double offset) { + return bin * cellSize + offset; + } + + /// Helper method to convert a 1D position to a cell ID + int Segmentation::positionToBin(double position, double cellSize, double offset) { + if (cellSize <= 1e-10) { + throw runtime_error("Invalid cell size: 0.0"); + } + return int(floor((position + 0.5 * cellSize - offset) / cellSize)); + } + + /// Helper method to convert a bin number to a 1D position given a vector of binBoundaries + double Segmentation::binToPosition(CellID bin, std::vector<double> const& cellBoundaries, double offset) { + return (cellBoundaries[bin+1] + cellBoundaries[bin])*0.5 + offset; + } + /// Helper method to convert a 1D position to a cell ID given a vector of binBoundaries + int Segmentation::positionToBin(double position, std::vector<double> const& cellBoundaries, double offset) { + + // include the lower edge to the segmentation, deal with numerical issues + if(fabs(position/cellBoundaries.front()-1.0) < 3e-12) return 0; + + // include the upper edge of the last bin to the segmentation, deal with numerical issues + if(fabs(position/cellBoundaries.back()-1.0) < 3e-12) return int(cellBoundaries.size()-2); + + // hits outside cannot be treated + if(position < cellBoundaries.front()) { + std::stringstream err; + err << std::setprecision(20) << std::scientific; + err << "Hit Position (" << position << ") is below the acceptance" + << " (" << cellBoundaries.front() << ") " + << "of the segmentation"; + throw std::runtime_error(err.str()); + } + if(position > cellBoundaries.back() ) { + std::stringstream err; + err << std::setprecision(20) << std::scientific; + err << "Hit Position (" << position << ") is above the acceptance" + << " (" << cellBoundaries.back() << ") " + << "of the segmentation"; + throw std::runtime_error(err.str()); + } + + + std::vector<double>::const_iterator bin = std::upper_bound(cellBoundaries.begin(), + cellBoundaries.end(), + position-offset); + + // need to reduce found bin by one, because upper_bound works that way, lower_bound would not help + return bin - cellBoundaries.begin() - 1 ; + + } + + std::vector<double> Segmentation::cellDimensions(const CellID&) const { + std::stringstream errorMessage; + errorMessage << __func__ << " is not implemented for " << _type; + throw std::logic_error(errorMessage.str()); + } + + } /* namespace DDSegmentation */ } /* namespace DD4hep */ diff --git a/cmake/DD4hepConfig.cmake.in b/cmake/DD4hepConfig.cmake.in index 5c5f127db..4b7697184 100644 --- a/cmake/DD4hepConfig.cmake.in +++ b/cmake/DD4hepConfig.cmake.in @@ -25,8 +25,8 @@ set ( DD4HEP_USE_CXX14 "@DD4HEP_USE_CXX14@" ) set ( Geant4_DIR "@Geant4_DIR@" ) set ( GEANT4_USE_CLHEP "@GEANT4_USE_CLHEP@" ) -set ( ROOT_DIR "@ROOTSYS@/cmake" ) -set ( ROOTSYS "@ROOTSYS@" ) +set ( ROOTSYS "$ENV{ROOTSYS}" ) +set ( ROOT_DIR "$ENV{ROOTSYS}/cmake" ) set ( ROOT_VERSION_MAJOR "@ROOT_VERSION_MAJOR@" ) include ( ${DD4hep_DIR}/cmake/DD4hep.cmake ) diff --git a/doc/release.notes b/doc/release.notes index 98737f819..12878e5b1 100644 --- a/doc/release.notes +++ b/doc/release.notes @@ -4,6 +4,22 @@ DD4hep ---- Release Notes ================================= +2016-06-24 M.Frank + Implement multiple segmentations. + Though one readout objects (associated to one DetElement) may only have on segmentation, + The MultiSegmentation type allows to have several sub-segmentations, which can be chosen + from. + Please see examples/ClientTests/*/MultiSegmentations + + At the same time allow the readout object to defined multiple collections through + the IDDescriptor. + Please see examples/ClientTests/*/MultiCollections + + The combined example can be found in + Please see examples/ClientTests/*/NestedSegmentations + + + 2016-05-03 M.Frank Green light is ON. You may use revision 2237 and above. diff --git a/examples/CLICSiD/compact/compact.xml b/examples/CLICSiD/compact/compact.xml index c6320d831..9395db50d 100644 --- a/examples/CLICSiD/compact/compact.xml +++ b/examples/CLICSiD/compact/compact.xml @@ -234,10 +234,10 @@ <vis name="HcalBarrelVis" alpha="0.1" r="1" g="1" b="0.1" showDaughters="true" visible="true"/> <vis name="HcalBarrelStavesVis" alpha="0.1" r="1" g="0" b="0.3" showDaughters="true" visible="true"/> <vis name="HcalBarrelLayerVis" alpha="0.5" r="1" g="0" b="0.5" showDaughters="true" visible="true"/> - <vis name="HcalBarrelSensorVis" alpha="1" r="1" g="1" b="0.7" showDaughters="true" visible="true"/> + <vis name="HcalBarrelSensorVis" alpha="1" r="1" g="1" b="0.7" showDaughters="true" visible="true"/> <vis name="HcalEndcapVis" alpha="0.1" r="1" g="1" b="0.1" showDaughters="false" visible="true"/> - <vis name="HcalEndcapLayerVis" alpha="1" r="1" g="0" b="0.5" showDaughters="false" visible="true"/> + <vis name="HcalEndcapLayerVis" alpha="1" r="1" g="0" b="0.5" showDaughters="false" visible="true"/> <vis name="SolenoidBarrelLayerVis" alpha="1" r="0" g="0.3" b="0.3" showDaughters="false" visible="true"/> <vis name="SolenoidCoilEndsVis" alpha="1" r="0" g="0.9" b="0.9" showDaughters="false" visible="true"/> diff --git a/examples/ClientTests/CMakeLists.txt b/examples/ClientTests/CMakeLists.txt index 3527c2263..e5a054df4 100644 --- a/examples/ClientTests/CMakeLists.txt +++ b/examples/ClientTests/CMakeLists.txt @@ -31,7 +31,8 @@ dd4hep_configure_scripts( ClientTests DEFAULT_SETUP WITH_TESTS) # # foreach (test Assemblies BoxTrafos IronCylinder LheD_tracker MagnetFields MaterialTester - MiniTel SectorBarrelCalorimeter SiliconBlock NestedSimple NestedDetectors ) + MiniTel SectorBarrelCalorimeter SiliconBlock NestedSimple NestedDetectors + MultiCollections MultiSegmentations ) # Test format conversions foreach( type lcdd gdml vis ) dd4hep_add_test_reg( ClientTests_${test}_converter_${type} @@ -48,7 +49,8 @@ endforeach() # IronCylinder has no segmentation! # MaterialTester no geometry # SectorBarrelCalorimeter is bad -foreach (test Assemblies BoxTrafos LheD_tracker MagnetFields MiniTel SiliconBlock NestedSimple NestedDetectors ) +foreach (test Assemblies BoxTrafos LheD_tracker MagnetFields MiniTel SiliconBlock + NestedSimple NestedDetectors MultiCollections ) # # Test material scans in [origine to 10 meters in y] dd4hep_add_test_reg( ClientTests_${test}_material_scan @@ -69,7 +71,7 @@ endforeach() # # # -foreach (test BoxTrafos IronCylinder MiniTel SiliconBlock NestedSimple ) +foreach (test BoxTrafos IronCylinder MiniTel SiliconBlock NestedSimple MultiCollections ) # # ROOT Geometry checks dd4hep_add_test_reg( ClientTests_${test}_check_geometry @@ -101,4 +103,14 @@ if (DD4HEP_USE_GEANT4) REGEX_PASS NONE REGEX_FAIL "Exception;EXCEPTION;ERROR;Error" ) endforeach(script) + # Geant4 full simulation checks of multi-collection/segmentation detectors + foreach(script MultiCollections MultiSegmentations MultiSegmentCollections ) + dd4hep_add_test_reg( ClientTests_sim_${script} + COMMAND "${CMAKE_INSTALL_PREFIX}/bin/run_test_ClientTests.sh" + EXEC_ARGS python ${CMAKE_CURRENT_SOURCE_DIR}/scripts/MultiCollections.py + -compact ${CMAKE_CURRENT_SOURCE_DIR}/compact/${script}.xml -batch + REQUIRES DDG4 Geant4 + REGEX_PASS NONE + REGEX_FAIL "Exception;EXCEPTION;ERROR;Error" ) + endforeach(script) endif(DD4HEP_USE_GEANT4) diff --git a/examples/ClientTests/compact/MultiCollections.xml b/examples/ClientTests/compact/MultiCollections.xml new file mode 100644 index 000000000..1b18f5cb5 --- /dev/null +++ b/examples/ClientTests/compact/MultiCollections.xml @@ -0,0 +1,64 @@ +<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0" + xmlns:xs="http://www.w3.org/2001/XMLSchema" + xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> + + <info name="Nested_Detectors_test" + title="Test for nested detector descriptions" + author="Markus Frank" + url="None" + status="development" + version="$Id: compact.xml 1374 2014-11-05 10:49:55Z markus.frank@cern.ch $"> + <comment>Test for nested detector descriptions</comment> + </info> + + <includes> + <gdmlFile ref="${DDDetectors_dir}/elements.xml"/> + <gdmlFile ref="${DDDetectors_dir}/materials.xml"/> + </includes> + + <define> + <constant name="world_side" value="30000*mm"/> + <constant name="world_x" value="world_side"/> + <constant name="world_y" value="world_side"/> + <constant name="world_z" value="world_side"/> + <constant name="DDDetectors_dir" value="${DD4hepINSTALL}/DDDetectors/compact" type="string"/>; + <constant name="SiD_dir" value="${DDDetectors_dir}/SiD" type="string"/>; + </define> + + <display> + <vis name="InvisibleNoDaughters" showDaughters="false" visible="false"/> + <vis name="InvisibleWithDaughters" showDaughters="true" visible="false"/> + <vis name="BlueVis" alpha="1" r="0.0" g="0.0" b="1.0" showDaughters="true" visible="true"/> + </display> + + <!-- ================================================================== --> + <!-- Simple calorimeter with a couple of layers --> + <!-- ================================================================== --> + <detectors> + <detector id="13" name="TestCal" reflect="true" type="DD4hep_CylindricalEndcapCalorimeter" readout="TestCalHits" vis="BlueVis"> + <comment>Test Calorimeter</comment> + <dimensions inner_r = "0.1*m" outer_r="2*m" inner_z = "2*m"/> + <layer repeat="15" > + <slice material = "Silicon" thickness = "0.032*cm" sensitive = "yes" /> + <slice material = "Copper" thickness = "0.005*cm" /> + <slice material = "Air" thickness = "0.033*cm" /> + </layer> + </detector> + </detectors> + + <!-- ================================================================== --> + <!-- Associated readout structures for the calorimeter --> + <!-- ================================================================== --> + <readouts> + <readout name="TestCalHits"> + <segmentation type="CartesianGridXY" grid_size_x="0.1*cm" grid_size_y="0.1*cm" /> + <hits_collections> + <hits_collection name="TestCallInnerLayerHits" key="layer" key_value="0x1"/> + <hits_collection name="TestCallMiddleLayerHits" key="layer" key_min="2" key_max="5"/> + <hits_collection name="TestCallOuterLayerHits" key="layer" key_min="0x6" key_max="0xFF"/> + </hits_collections> + <id>system:8,barrel:3,layer:8,slice:8,x:32:-16,y:-16</id> + </readout> + </readouts> + +</lccdd> diff --git a/examples/ClientTests/compact/MultiSegmentCollections.xml b/examples/ClientTests/compact/MultiSegmentCollections.xml new file mode 100644 index 000000000..599c8c598 --- /dev/null +++ b/examples/ClientTests/compact/MultiSegmentCollections.xml @@ -0,0 +1,67 @@ +<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0" + xmlns:xs="http://www.w3.org/2001/XMLSchema" + xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> + + <info name="Nested_Detectors_test" + title="Test for nested detector descriptions" + author="Markus Frank" + url="None" + status="development" + version="$Id: compact.xml 1374 2014-11-05 10:49:55Z markus.frank@cern.ch $"> + <comment>Test for nested detector descriptions</comment> + </info> + + <includes> + <gdmlFile ref="${DDDetectors_dir}/elements.xml"/> + <gdmlFile ref="${DDDetectors_dir}/materials.xml"/> + </includes> + + <define> + <constant name="world_side" value="30000*mm"/> + <constant name="world_x" value="world_side"/> + <constant name="world_y" value="world_side"/> + <constant name="world_z" value="world_side"/> + <constant name="DDDetectors_dir" value="${DD4hepINSTALL}/DDDetectors/compact" type="string"/>; + <constant name="SiD_dir" value="${DDDetectors_dir}/SiD" type="string"/>; + </define> + + <display> + <vis name="InvisibleNoDaughters" showDaughters="false" visible="false"/> + <vis name="InvisibleWithDaughters" showDaughters="true" visible="false"/> + <vis name="BlueVis" alpha="1" r="0.0" g="0.0" b="1.0" showDaughters="true" visible="true"/> + </display> + + <!-- ================================================================== --> + <!-- Simple calorimeter with a couple of layers --> + <!-- ================================================================== --> + <detectors> + <detector id="13" name="TestCal" reflect="true" type="DD4hep_CylindricalEndcapCalorimeter" readout="TestCalHits" vis="BlueVis"> + <comment>Test Calorimeter</comment> + <dimensions inner_r = "0.1*m" outer_r="2*m" inner_z = "2*m"/> + <layer repeat="15" > + <slice material = "Silicon" thickness = "0.032*cm" sensitive = "yes" /> + <slice material = "Copper" thickness = "0.005*cm" /> + <slice material = "Air" thickness = "0.033*cm" /> + </layer> + </detector> + </detectors> + + <!-- ================================================================== --> + <!-- Associated readout structures for the calorimeter --> + <!-- ================================================================== --> + <readouts> + <readout name="TestCalHits"> + <segmentation type="MultiSegmentation" key="layer"> + <segmentation name="Layer1grid" type="CartesianGridXY" key_min="0x1" key_max="4" grid_size_x="0.1*cm" grid_size_y="0.1*cm" /> + <segmentation name="Layer2grid" type="CartesianGridXY" key_value="5" grid_size_x="0.2*cm" grid_size_y="0.2*cm" /> + <segmentation name="Layer3grid" type="CartesianGridXY" key_min="0x6" key_max="0xFF" grid_size_x="0.3*cm" grid_size_y="0.3*cm" /> + </segmentation> + <hits_collections> + <hits_collection name="TestCallInnerLayerHits" key="layer" key_value="0x1"/> + <hits_collection name="TestCallMiddleLayerHits" key="layer" key_min="2" key_max="5"/> + <hits_collection name="TestCallOuterLayerHits" key="layer" key_min="0x6" key_max="0xFF"/> + </hits_collections> + <id>system:8,barrel:3,layer:8,slice:8,x:32:-16,y:-16</id> + </readout> + </readouts> +</lccdd> diff --git a/examples/ClientTests/compact/MultiSegmentations.xml b/examples/ClientTests/compact/MultiSegmentations.xml new file mode 100644 index 000000000..2e086697f --- /dev/null +++ b/examples/ClientTests/compact/MultiSegmentations.xml @@ -0,0 +1,62 @@ +<lccdd xmlns:compact="http://www.lcsim.org/schemas/compact/1.0" + xmlns:xs="http://www.w3.org/2001/XMLSchema" + xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> + + <info name="Nested_Detectors_test" + title="Test for nested detector descriptions" + author="Markus Frank" + url="None" + status="development" + version="$Id: compact.xml 1374 2014-11-05 10:49:55Z markus.frank@cern.ch $"> + <comment>Test for nested detector descriptions</comment> + </info> + + <includes> + <gdmlFile ref="${DDDetectors_dir}/elements.xml"/> + <gdmlFile ref="${DDDetectors_dir}/materials.xml"/> + </includes> + + <define> + <constant name="world_side" value="30000*mm"/> + <constant name="world_x" value="world_side"/> + <constant name="world_y" value="world_side"/> + <constant name="world_z" value="world_side"/> + <constant name="DDDetectors_dir" value="${DD4hepINSTALL}/DDDetectors/compact" type="string"/>; + <constant name="SiD_dir" value="${DDDetectors_dir}/SiD" type="string"/>; + </define> + + <display> + <vis name="InvisibleNoDaughters" showDaughters="false" visible="false"/> + <vis name="InvisibleWithDaughters" showDaughters="true" visible="false"/> + <vis name="BlueVis" alpha="1" r="0.0" g="0.0" b="1.0" showDaughters="true" visible="true"/> + </display> + + <!-- ================================================================== --> + <!-- Simple calorimeter with a couple of layers --> + <!-- ================================================================== --> + <detectors> + <detector id="13" name="TestCal" reflect="true" type="DD4hep_CylindricalEndcapCalorimeter" readout="TestCalHits" vis="BlueVis"> + <comment>Test Calorimeter</comment> + <dimensions inner_r = "0.1*m" outer_r="2*m" inner_z = "2*m"/> + <layer repeat="15" > + <slice material = "Silicon" thickness = "0.032*cm" sensitive = "yes" /> + <slice material = "Copper" thickness = "0.005*cm" /> + <slice material = "Air" thickness = "0.033*cm" /> + </layer> + </detector> + </detectors> + + <!-- ================================================================== --> + <!-- Associated readout structures for the calorimeter --> + <!-- ================================================================== --> + <readouts> + <readout name="TestCalHits"> + <segmentation type="MultiSegmentation" key="layer"> + <segmentation name="Layer1grid" type="CartesianGridXY" key_min="0x1" key_max="4" grid_size_x="0.1*cm" grid_size_y="0.1*cm" /> + <segmentation name="Layer2grid" type="CartesianGridXY" key_value="5" grid_size_x="0.2*cm" grid_size_y="0.2*cm" /> + <segmentation name="Layer3grid" type="CartesianGridXY" key_min="0x6" key_max="0xFF" grid_size_x="0.3*cm" grid_size_y="0.3*cm" /> + </segmentation> + <id>system:8,barrel:3,layer:8,slice:8,x:32:-16,y:-16</id> + </readout> + </readouts> +</lccdd> diff --git a/examples/ClientTests/scripts/MultiCollections.py b/examples/ClientTests/scripts/MultiCollections.py new file mode 100644 index 000000000..03b9d0fc7 --- /dev/null +++ b/examples/ClientTests/scripts/MultiCollections.py @@ -0,0 +1,62 @@ +import os, sys, time, DDG4 +from DDG4 import OutputLevel as Output +from SystemOfUnits import * +# +# +""" + + DD4hep example setup using the python configuration + + \author M.Frank + \version 1.0 + +""" +def run(): + batch = False + kernel = DDG4.Kernel() + install_dir = os.environ['DD4hepINSTALL'] + geometry = "file:"+install_dir+"/examples/ClientTests/compact/MultiCollections.xml" + kernel.setOutputLevel('Geant4Converter',Output.DEBUG) + kernel.setOutputLevel('Gun',Output.INFO) + for i in xrange(len(sys.argv)): + if sys.argv[i]=='-compact': + geometry = sys.argv[i+1] + elif sys.argv[i]=='-input': + geometry = sys.argv[i+1] + elif sys.argv[i]=='-batch': + batch = True + elif sys.argv[i]=='batch': + batch = True + + kernel.loadGeometry(geometry) + geant4 = DDG4.Geant4(kernel) + geant4.printDetectors() + geant4.setupCshUI() + if batch: kernel.UI = '' + + # Configure field + field = geant4.setupTrackingField(prt=True) + # Configure I/O + evt_root = geant4.setupROOTOutput('RootOutput','Multi_coll_'+time.strftime('%Y-%m-%d_%H-%M'),mc_truth=True) + # Setup particle gun + geant4.setupGun("Gun",particle='pi-',energy=10*GeV,multiplicity=1) + + # Now the test calorimeter with multiple collections + seq,act = geant4.setupCalorimeter('TestCal') + + # And handle the simulation particles. + part = DDG4.GeneratorAction(kernel,"Geant4ParticleHandler/ParticleHandler") + kernel.generatorAction().adopt(part) + part.MinimalKineticEnergy = 1*MeV + part.enableUI() + + # Now build the physics list: + phys = kernel.physicsList() + phys.extends = 'QGSP_BERT' + phys.enableUI() + phys.dump() + # and run + geant4.execute() + +if __name__ == "__main__": + run() -- GitLab