From 7ffff51e18750299916f3afd864a6dfddc2d65b1 Mon Sep 17 00:00:00 2001 From: Markus Frank <Markus.Frank@cern.ch> Date: Fri, 8 Mar 2019 11:10:45 +0100 Subject: [PATCH] Towards a first Geant4 simulation of optical photons -- now have to wait for TGeo add-ons --- DDCore/include/DD4hep/OpticalSurfaceManager.h | 18 +- DDCore/include/DD4hep/OpticalSurfaces.h | 21 + DDCore/src/OpticalSurfaceManager.cpp | 21 +- DDCore/src/plugins/Compact2Objects.cpp | 61 ++- DDCore/src/plugins/StandardPlugins.cpp | 97 ++++- DDG4/include/DDG4/Geant4Converter.h | 2 + DDG4/include/DDG4/Geant4PhysicsList.h | 19 +- DDG4/plugins/Geant4CerenkovPhysics.cpp | 108 +++++ .../Geant4DetectorGeometryConstruction.cpp | 4 + DDG4/plugins/Geant4OpticalPhotonPhysics.cpp | 104 +++++ DDG4/plugins/Geant4ScintillationPhysics.cpp | 119 ++++++ DDG4/src/Geant4Converter.cpp | 46 ++- DDG4/src/Geant4PhysicsList.cpp | 151 +++++-- examples/OpticalSurfaces/CMakeLists.txt | 16 +- examples/OpticalSurfaces/compact/OpNovice.xml | 375 ++++++++++++++++++ .../compact/ReadMaterialProperties.xml | 233 +++++++++++ .../compact/ReadOpticalSurfaces.xml | 21 +- examples/OpticalSurfaces/scripts/OpNovice.py | 116 ++++++ examples/OpticalSurfaces/src/OpNovice_geo.cpp | 58 +++ 19 files changed, 1487 insertions(+), 103 deletions(-) create mode 100644 DDG4/plugins/Geant4CerenkovPhysics.cpp create mode 100644 DDG4/plugins/Geant4OpticalPhotonPhysics.cpp create mode 100644 DDG4/plugins/Geant4ScintillationPhysics.cpp create mode 100644 examples/OpticalSurfaces/compact/OpNovice.xml create mode 100644 examples/OpticalSurfaces/compact/ReadMaterialProperties.xml create mode 100644 examples/OpticalSurfaces/scripts/OpNovice.py create mode 100644 examples/OpticalSurfaces/src/OpNovice_geo.cpp diff --git a/DDCore/include/DD4hep/OpticalSurfaceManager.h b/DDCore/include/DD4hep/OpticalSurfaceManager.h index ccefe96de..a5029c112 100644 --- a/DDCore/include/DD4hep/OpticalSurfaceManager.h +++ b/DDCore/include/DD4hep/OpticalSurfaceManager.h @@ -59,12 +59,18 @@ namespace dd4hep { static OpticalSurfaceManager getOpticalSurfaceManager(Detector& description); #if ROOT_VERSION_CODE > ROOT_VERSION(6,16,0) - /// Access skin surface by its identifier - SkinSurface getSkinSurface(DetElement de, const std::string& nam) const; - /// Access border surface by its identifier - BorderSurface getBorderSurface(DetElement de, const std::string& nam) const; - /// Access optical surface data by its identifier - OpticalSurface getOpticalSurface(DetElement de, const std::string& nam) const; + /// Access skin surface by its full name + SkinSurface skinSurface(const std::string& full_name) const; + /// Access skin surface by its identifier tuple (DetElement, name) + SkinSurface skinSurface(DetElement de, const std::string& nam) const; + /// Access border surface by its full name + BorderSurface borderSurface(const std::string& full_name) const; + /// Access border surface by its identifier tuple (DetElement, name) + BorderSurface borderSurface(DetElement de, const std::string& nam) const; + /// Access optical surface data by its full name + OpticalSurface opticalSurface(const std::string& full_name) const; + /// Access optical surface data by its identifier tuple (DetElement, name) + OpticalSurface opticalSurface(DetElement de, const std::string& nam) const; /// Add skin surface to manager void addSkinSurface(DetElement de, SkinSurface surf) const; /// Add border surface to manager diff --git a/DDCore/include/DD4hep/OpticalSurfaces.h b/DDCore/include/DD4hep/OpticalSurfaces.h index 77487a67a..0cf4c7540 100644 --- a/DDCore/include/DD4hep/OpticalSurfaces.h +++ b/DDCore/include/DD4hep/OpticalSurfaces.h @@ -73,6 +73,27 @@ namespace dd4hep { /// Assignment operator OpticalSurface& operator=(const OpticalSurface& m) = default; + + /// Convenience function forwarding to TGeoOpticalSurface + static EType stringToType(const std::string& type) { + return Object::StringToType(type.c_str()); + } + /// Convenience function forwarding to TGeoOpticalSurface + static std::string typeToString(EType type) { + return Object::TypeToString(type); + } + static EModel stringToModel(const std::string& model) { + return Object::StringToModel(model.c_str()); + } + static std::string modelToString(EModel model) { + return Object::ModelToString(model); + } + static EFinish stringToFinish(const std::string& finish) { + return Object::StringToFinish(finish.c_str()); + } + static std::string finishToString(EFinish finish) { + return Object::FinishToString(finish); + } }; /// Class to support the handling of optical surfaces. diff --git a/DDCore/src/OpticalSurfaceManager.cpp b/DDCore/src/OpticalSurfaceManager.cpp index 43a97ac49..51259154e 100644 --- a/DDCore/src/OpticalSurfaceManager.cpp +++ b/DDCore/src/OpticalSurfaceManager.cpp @@ -33,7 +33,7 @@ OpticalSurfaceManager OpticalSurfaceManager::getOpticalSurfaceManager(Detector& #if ROOT_VERSION_CODE > ROOT_VERSION(6,16,0) /// Access skin surface by its identifier -SkinSurface OpticalSurfaceManager::getSkinSurface(DetElement de, const string& nam) const { +SkinSurface OpticalSurfaceManager::skinSurface(DetElement de, const string& nam) const { if ( de.isValid() ) { Object* o = access(); string n = de.path() + '#' + nam; @@ -48,8 +48,13 @@ SkinSurface OpticalSurfaceManager::getSkinSurface(DetElement de, const string& return SkinSurface(); } +/// Access skin surface by its full name +SkinSurface OpticalSurfaceManager::skinSurface(const string& full_nam) const { + return access()->detector.manager().GetSkinSurface(full_nam.c_str()); +} + /// Access border surface by its identifier -BorderSurface OpticalSurfaceManager::getBorderSurface(DetElement de, const string& nam) const { +BorderSurface OpticalSurfaceManager::borderSurface(DetElement de, const string& nam) const { if ( de.isValid() ) { Object* o = access(); string n = de.path() + '#' + nam; @@ -64,8 +69,13 @@ BorderSurface OpticalSurfaceManager::getBorderSurface(DetElement de, const stri return BorderSurface(); } +/// Access border surface by its full name +BorderSurface OpticalSurfaceManager::borderSurface(const string& full_nam) const { + return access()->detector.manager().GetBorderSurface(full_nam.c_str()); +} + /// Access optical surface data by its identifier -OpticalSurface OpticalSurfaceManager::getOpticalSurface(DetElement de, const string& nam) const { +OpticalSurface OpticalSurfaceManager::opticalSurface(DetElement de, const string& nam) const { if ( de.isValid() ) { Object* o = access(); string n = de.path() + '#' + nam; @@ -80,6 +90,11 @@ OpticalSurface OpticalSurfaceManager::getOpticalSurface(DetElement de, const str return OpticalSurface(); } +/// Access optical surface data by its identifier +OpticalSurface OpticalSurfaceManager::opticalSurface(const string& full_nam) const { + return access()->detector.manager().GetOpticalSurface(full_nam.c_str()); +} + /// Add skin surface to manager void OpticalSurfaceManager::addSkinSurface(DetElement de, SkinSurface surf) const { if ( access()->skinSurfaces.insert(make_pair(make_pair(de,surf->GetName()), surf)).second ) diff --git a/DDCore/src/plugins/Compact2Objects.cpp b/DDCore/src/plugins/Compact2Objects.cpp index e2ff3f86c..2996f4210 100644 --- a/DDCore/src/plugins/Compact2Objects.cpp +++ b/DDCore/src/plugins/Compact2Objects.cpp @@ -33,6 +33,7 @@ // Root/TGeo include files #include "TGeoManager.h" #include "TGeoMaterial.h" +#include "TGDMLMatrix.h" // C/C++ include files #include <climits> @@ -465,8 +466,12 @@ template <> void Converter<Material>::operator()(xml_h e) const { string ref = p.attr<string>(_U(ref)); TGDMLMatrix* m = mgr.GetGDMLMatrix(ref.c_str()); if ( m ) { + string prop_nam = p.attr<string>(_U(name)); //TODO: mat->AddProperty(p.attr<string>(_U(name)).c_str(), m); - mat->AddProperty(p.attr<string>(_U(name)).c_str(), ref.c_str()); + mat->AddProperty(prop_nam.c_str(), ref.c_str()); + printout(s_debug.materials ? ALWAYS : DEBUG, "Compact", + "++ material %-16s add property: %s -> %s.", + mat->GetName(), prop_nam.c_str(), ref.c_str()); continue; } // ERROR @@ -618,27 +623,50 @@ template <> void Converter<Atom>::operator()(xml_h e) const { * */ template <> void Converter<OpticalSurface>::operator()(xml_h element) const { - xml_attr_t attr; xml_elt_t e = element; + // Defaults from Geant4 OpticalSurface::EModel model = OpticalSurface::Model::kMglisur; OpticalSurface::EFinish finish = OpticalSurface::Finish::kFpolished; OpticalSurface::EType type = OpticalSurface::Type::kTdielectric_metal; - Double_t value = 0; - if ( (attr=e.attr<xml_attr_t>(_U(type))) ) type = (OpticalSurface::EType)e.attr<int>(attr); - if ( (attr=e.attr<xml_attr_t>(_U(model))) ) model = (OpticalSurface::EModel)e.attr<int>(attr); - if ( (attr=e.attr<xml_attr_t>(_U(finish))) ) finish = (OpticalSurface::EFinish)e.attr<int>(attr); - if ( (attr=e.attr<xml_attr_t>(_U(value))) ) value = e.attr<double>(attr); + Double_t value = 1.0; + if ( e.hasAttr(_U(type)) ) type = OpticalSurface::stringToType(e.attr<string>(_U(type))); + if ( e.hasAttr(_U(model)) ) model = OpticalSurface::stringToModel(e.attr<string>(_U(model))); + if ( e.hasAttr(_U(finish)) ) finish = OpticalSurface::stringToFinish(e.attr<string>(_U(finish))); + if ( e.hasAttr(_U(value)) ) value = e.attr<double>(_U(value)); OpticalSurface surf(description, e.attr<string>(_U(name)), model, finish, type, value); if ( s_debug.surface ) { printout(ALWAYS,"Compact","+++ Reading optical surface %s Typ:%d model:%d finish:%d value:%f", e.attr<string>(_U(name)).c_str(), int(type), int(model), int(finish), value); } for (xml_coll_t props(e, _U(property)); props; ++props) { - surf->AddProperty(props.attr<string>(_U(name)).c_str(), props.attr<string>(_U(ref)).c_str()); - if ( s_debug.surface ) { - printout(ALWAYS,"Compact","+++ \t\t Property: %s -> %s", - props.attr<string>(_U(name)).c_str(), props.attr<string>(_U(ref)).c_str()); + if ( props.hasAttr(_U(ref)) ) { + surf->AddProperty(props.attr<string>(_U(name)).c_str(), props.attr<string>(_U(ref)).c_str()); + if ( s_debug.surface ) { + printout(ALWAYS,"Compact","+++ \t\t Property: %s -> %s", + props.attr<string>(_U(name)).c_str(), props.attr<string>(_U(ref)).c_str()); + } + continue; + } + size_t cols = props.attr<long>(_U(coldim)); + string nam = props.attr<string>(_U(name)); + stringstream str(props.attr<string>(_U(values))), str_nam; + string val; + vector<double> values; + while ( !str.eof() ) { + val = ""; + str >> val; + if ( val.empty() && !str.good() ) break; + values.push_back(_toDouble(val)); } + /// Create table and register table + TGDMLMatrix* table = new TGDMLMatrix("",values.size()/cols, cols); + str_nam << nam << "__" << (void*)table; + table->SetName(str_nam.str().c_str()); + table->SetTitle(nam.c_str()); + for (size_t i=0, n=values.size(); i<n; ++i) + table->Set(i/cols, i%cols, values[i]); + surf->AddProperty(nam.c_str(), table->GetName()); + description.manager().AddGDMLMatrix(table); } } @@ -659,13 +687,14 @@ template <> void Converter<PropertyTable>::operator()(xml_h e) const { } values.reserve(1024); while ( !str.eof() ) { + val = ""; str >> val; - if ( !str.good() ) break; + if ( val.empty() && !str.good() ) break; + values.push_back(_toDouble(val)); if ( s_debug.matrix ) { cout << " state:" << (str.good() ? "OK " : "BAD") << " '" << val << "'"; + if ( 0 == (values.size()%cols) ) cout << endl; } - values.push_back(_toDouble(val)); - if ( 0 == (values.size()%cols) ) cout << endl; } if ( s_debug.matrix ) { cout << endl; @@ -1370,14 +1399,14 @@ template <> void Converter<Compact>::operator()(xml_h element) const { if (element.hasChild(_U(info))) (Converter<Header>(description))(xml_h(compact.child(_U(info)))); - xml_coll_t(compact, _U(materials)).for_each(_U(element), Converter<Atom>(description)); - xml_coll_t(compact, _U(materials)).for_each(_U(material), Converter<Material>(description)); xml_coll_t(compact, _U(properties)).for_each(_U(attributes), Converter<Property>(description)); #if ROOT_VERSION_CODE > ROOT_VERSION(6,16,0) /// These two must be parsed early, because they are needed by the detector constructors xml_coll_t(compact, _U(properties)).for_each(_U(matrix), Converter<PropertyTable>(description)); xml_coll_t(compact, _U(surfaces)).for_each(_U(opticalsurface), Converter<OpticalSurface>(description)); #endif + xml_coll_t(compact, _U(materials)).for_each(_U(element), Converter<Atom>(description)); + xml_coll_t(compact, _U(materials)).for_each(_U(material), Converter<Material>(description)); xml_coll_t(compact, _U(display)).for_each(_U(include), Converter<DetElementInclude>(description)); xml_coll_t(compact, _U(display)).for_each(_U(vis), Converter<VisAttr>(description)); diff --git a/DDCore/src/plugins/StandardPlugins.cpp b/DDCore/src/plugins/StandardPlugins.cpp index 4583514c3..8495b3701 100644 --- a/DDCore/src/plugins/StandardPlugins.cpp +++ b/DDCore/src/plugins/StandardPlugins.cpp @@ -306,10 +306,11 @@ DECLARE_APPLY(DD4hep_Dump_GDMLTables,root_dump_gdml_tables) */ static long root_dump_optical_surfaces(Detector& description, int /* argc */, char** /* argv */) { size_t num_surfaces = 0; + printout(ALWAYS,"",""); #if ROOT_VERSION_CODE > ROOT_VERSION(6,16,0) const TObjArray* c = description.manager().GetListOfOpticalSurfaces(); TObjArrayIter arr(c); - printout(INFO,"Dump_OpticalSurfaces","+++ Dumping known Optical Surfaces from TGeoManager."); + printout(ALWAYS,"Dump_OpticalSurfaces","+++ Dumping known Optical Surfaces from TGeoManager."); for(TObject* i = arr.Next(); i; i=arr.Next()) { TGeoOpticalSurface* m = (TGeoOpticalSurface*)i; ++num_surfaces; @@ -318,12 +319,72 @@ static long root_dump_optical_surfaces(Detector& description, int /* argc */, ch #else description.world().isValid(); #endif - printout(INFO,"Dump_OpticalSurfaces", + printout(ALWAYS,"Dump_OpticalSurfaces", "+++ Successfully dumped %ld Optical surfaces.",num_surfaces); return 1; } DECLARE_APPLY(DD4hep_Dump_OpticalSurfaces,root_dump_optical_surfaces) +/// Basic entry point to dump all known skin surfaces +/** + * + * Factory: DD4hep_Dump_SkinSurfaces + * + * \author M.Frank + * \version 1.0 + * \date 01/04/2014 + */ +static long root_dump_skin_surfaces(Detector& description, int /* argc */, char** /* argv */) { + size_t num_surfaces = 0; + printout(ALWAYS,"",""); +#if ROOT_VERSION_CODE > ROOT_VERSION(6,16,0) + const TObjArray* c = description.manager().GetListOfSkinSurfaces(); + TObjArrayIter arr(c); + printout(ALWAYS,"Dump_SkinSurfaces","+++ Dumping known Skin Surfaces from TGeoManager."); + for(TObject* i = arr.Next(); i; i=arr.Next()) { + TGeoSkinSurface* m = (TGeoSkinSurface*)i; + ++num_surfaces; + m->Print(); + } +#else + description.world().isValid(); +#endif + printout(ALWAYS,"Dump_SkinSurfaces", + "+++ Successfully dumped %ld Skin surfaces.",num_surfaces); + return 1; +} +DECLARE_APPLY(DD4hep_Dump_SkinSurfaces,root_dump_skin_surfaces) + +/// Basic entry point to dump all known border surfaces +/** + * + * Factory: DD4hep_Dump_BorderSurfaces + * + * \author M.Frank + * \version 1.0 + * \date 01/04/2014 + */ +static long root_dump_border_surfaces(Detector& description, int /* argc */, char** /* argv */) { + size_t num_surfaces = 0; + printout(ALWAYS,"",""); +#if ROOT_VERSION_CODE > ROOT_VERSION(6,16,0) + const TObjArray* c = description.manager().GetListOfBorderSurfaces(); + TObjArrayIter arr(c); + printout(ALWAYS,"Dump_BorderSurfaces","+++ Dumping known Border Surfaces from TGeoManager."); + for(TObject* i = arr.Next(); i; i=arr.Next()) { + TGeoBorderSurface* m = (TGeoBorderSurface*)i; + ++num_surfaces; + m->Print(); + } +#else + description.world().isValid(); +#endif + printout(ALWAYS,"Dump_BorderSurfaces", + "+++ Successfully dumped %ld Border surfaces.",num_surfaces); + return 1; +} +DECLARE_APPLY(DD4hep_Dump_BorderSurfaces,root_dump_border_surfaces) + /// Basic entry point to dump the ROOT TGeoElementTable object /** * Dump the elment table to stdout or file. @@ -442,7 +503,8 @@ DECLARE_APPLY(DD4hep_ElementTable,root_elements) static long root_materials(Detector& description, int argc, char** argv) { struct MaterialPrint { typedef xml::Element elt_h; - MaterialPrint() = default; + Detector& description; + MaterialPrint(Detector& desc) : description(desc) {} virtual ~MaterialPrint() = default; virtual elt_h print(TGeoMaterial* mat) { ::printf("%-8s %-32s Aeff=%7.3f Zeff=%7.4f rho=%8.3f [g/mole] radlen=%8.3g [cm] intlen=%8.3g [cm] index=%3d\n", @@ -455,6 +517,14 @@ static long root_materials(Detector& description, int argc, char** argv) { ::printf(" %-6s Fraction: %7.3f Z=%3d A=%6.2f N=%3d Neff=%6.2f\n", elt->GetName(), frac, elt->Z(), elt->A(), elt->N(), elt->Neff()); } + virtual void printProperty(elt_h, TNamed* prop, TGDMLMatrix* matrix) { + if ( matrix ) + ::printf(" Property: %-20s [%ld x %ld] --> %s\n", + prop->GetName(), long(matrix->GetRows()), long(matrix->GetCols()), prop->GetTitle()); + else + ::printf(" Property: %-20s [ERROR: NO TABLE!] --> %s\n", + prop->GetName(), prop->GetTitle()); + } virtual void operator()(TGeoMaterial* mat) { Double_t* mix = mat->IsMixture() ? ((TGeoMixture*)mat)->GetWmixt() : 0; elt_h mh = print(mat); @@ -462,11 +532,15 @@ static long root_materials(Detector& description, int argc, char** argv) { TGeoElement* elt = mat->GetElement(i); print(mh, elt, mix ? mix[i] : 1); } + TListIter mat_iter(&mat->GetProperties()); + for(TObject* i = mat_iter.Next(); i; i=mat_iter.Next()) { + printProperty(mh, (TNamed*)i, description.manager().GetGDMLMatrix(i->GetTitle())); + } } }; struct MaterialPrintXML : public MaterialPrint { elt_h root; - MaterialPrintXML(elt_h r) : root(r) {} + MaterialPrintXML(elt_h elt, Detector& desc) : MaterialPrint(desc), root(elt) {} virtual ~MaterialPrintXML() {} virtual elt_h print(TGeoMaterial* mat) { elt_h elt = root.addChild(_U(material)); @@ -488,13 +562,19 @@ static long root_materials(Detector& description, int argc, char** argv) { elt.setAttr(_U(n),frac); elt.setAttr(_U(ref),element->GetName()); } + virtual void printProperty(elt_h mat, TNamed* prop, TGDMLMatrix* /* matrix */) { + elt_h elt = mat.addChild(_U(property)); + elt.setAttr(_U(name),prop->GetName()); + elt.setAttr(_U(ref), prop->GetTitle()); + } }; - string type = "text", output = ""; + string type = "text", output = "", name = ""; for(int i=0; i<argc; ++i) { if ( argv[i][0] == '-' ) { char c = ::tolower(argv[i][1]); if ( c == 't' && i+1<argc ) type = argv[++i]; + else if ( c == 'n' && i+1<argc ) name = argv[++i]; else if ( c == 'o' && i+1<argc ) output = argv[++i]; else { ::printf("DD4hep_MaterialTable -opt [-opt] \n" @@ -527,14 +607,15 @@ static long root_materials(Detector& description, int argc, char** argv) { element = doc.root(); } dd4hep_ptr<MaterialPrint> printer(element - ? new MaterialPrintXML(element) - : new MaterialPrint()); + ? new MaterialPrintXML(element, description) + : new MaterialPrint(description)); TObject* obj = 0; TList* mats = description.manager().GetListOfMaterials(); dd4hep_ptr<TIterator> iter(mats->MakeIterator()); while( (obj=iter->Next()) != 0 ) { TGeoMaterial* mat = (TGeoMaterial*)obj; - (*printer)(mat); + if ( name.empty() || name == mat->GetName() ) + (*printer)(mat); } if ( element ) { xml::DocumentHandler dH; diff --git a/DDG4/include/DDG4/Geant4Converter.h b/DDG4/include/DDG4/Geant4Converter.h index b3fd7ffa9..71478cff1 100644 --- a/DDG4/include/DDG4/Geant4Converter.h +++ b/DDG4/include/DDG4/Geant4Converter.h @@ -43,6 +43,8 @@ namespace dd4hep { bool debugPlacements = false; /// Property: Flag to debug regions during conversion mechanism bool debugRegions = false; + /// Property: Flag to debug surfaces during conversion mechanism + bool debugSurfaces = false; /// Property: Flag to dump all placements after the conversion procedure bool printPlacements = false; diff --git a/DDG4/include/DDG4/Geant4PhysicsList.h b/DDG4/include/DDG4/Geant4PhysicsList.h index 91aea5d47..6875fc102 100644 --- a/DDG4/include/DDG4/Geant4PhysicsList.h +++ b/DDG4/include/DDG4/Geant4PhysicsList.h @@ -53,7 +53,7 @@ namespace dd4hep { class Process { public: std::string name; - int ordAtRestDoIt, ordAlongSteptDoIt, ordPostStepDoIt; + int ordAtRestDoIt=-1, ordAlongSteptDoIt=-1, ordPostStepDoIt=-1; /// Default constructor Process(); /// Copy constructor @@ -101,6 +101,7 @@ namespace dd4hep { typedef std::vector<PhysicsConstructor> PhysicsConstructors; PhysicsProcesses m_processes; + PhysicsProcesses m_discreteProcesses; PhysicsConstructors m_physics; ParticleConstructors m_particles; ParticleConstructors m_particlegroups; @@ -126,6 +127,20 @@ namespace dd4hep { ParticleProcesses& processes(const std::string& part_name); /// Access processes for one particle type (CONST) const ParticleProcesses& processes(const std::string& part_name) const; + + /// Access all physics discrete processes + PhysicsProcesses& discreteProcesses() { + return m_discreteProcesses; + } + /// Access all physics discrete processes + const PhysicsProcesses& discreteProcesses() const { + return m_discreteProcesses; + } + /// Access discrete processes for one particle type + ParticleProcesses& discreteProcesses(const std::string& part_name); + /// Access discrete processes for one particle type (CONST) + const ParticleProcesses& discreteProcesses(const std::string& part_name) const; + /// Access all physics particles ParticleConstructors& particles() { return m_particles; @@ -158,6 +173,8 @@ namespace dd4hep { /// Add particle process by name with arguments void addParticleProcess(const std::string& part_name, const std::string& proc_name, int ordAtRestDoIt,int ordAlongSteptDoIt,int ordPostStepDoIt); + /// Add discrete particle process by name with arguments + void addDiscreteParticleProcess(const std::string& part_name, const std::string& proc_name); /// Add PhysicsConstructor by name /** This constructor is used for intrinsic Geant4 defined physics constructors. * Such physics constructors are only created by the factory and attached diff --git a/DDG4/plugins/Geant4CerenkovPhysics.cpp b/DDG4/plugins/Geant4CerenkovPhysics.cpp new file mode 100644 index 000000000..a2eedcf3f --- /dev/null +++ b/DDG4/plugins/Geant4CerenkovPhysics.cpp @@ -0,0 +1,108 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// 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 +// +//========================================================================== + +/** \addtogroup Geant4PhysicsConstructor + * + * @{ + * \package Geant4CerenkovPhysics + * \brief PhysicsConstructor for Cerenkov physics + + This plugin allows to enable Cerenkov physics in the physics list + + * @} + */ +#ifndef DDG4_GEANT4CERENKOVPHYSICS_H +#define DDG4_GEANT4CERENKOVPHYSICS_H 1 + +/// Framework include files +#include "DDG4/Geant4PhysicsList.h" + +/// Geant4 include files +#include "G4ParticleTableIterator.hh" +#include "G4ParticleDefinition.hh" +#include "G4ParticleTypes.hh" +#include "G4ParticleTable.hh" +#include "G4ProcessManager.hh" + +#include "G4Cerenkov.hh" + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace sim { + + /// Geant4 physics list action to enable Cerenkov physics + /** + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_SIMULATION + */ + class Geant4CerenkovPhysics : public Geant4PhysicsList { + public: + /// Default constructor + Geant4CerenkovPhysics() = delete; + /// Copy constructor + Geant4CerenkovPhysics(const Geant4CerenkovPhysics&) = delete; + /// Initializing constructor + Geant4CerenkovPhysics(Geant4Context* ctxt, const std::string& nam) + : Geant4PhysicsList(ctxt, nam) + { + declareProperty("MaxNumPhotonsPerStep", m_maxNumPhotonsPerStep = 0); + declareProperty("MaxBetaChangePerStep", m_maxBetaChangePerStep = 0.0); + declareProperty("TrackSecondariesFirst", m_trackSecondariesFirst = false); + declareProperty("StackPhotons", m_stackPhotons = true); + declareProperty("VerboseLevel", m_verbosity = 0); + } + /// Default destructor + virtual ~Geant4CerenkovPhysics() = default; + /// Callback to construct processes (uses the G4 particle table) + virtual void constructProcesses(G4VUserPhysicsList* physics_list) { + this->Geant4PhysicsList::constructProcesses(physics_list); + info("+++ Constructing: maxNumPhotonsPerStep:%d maxBeta:%f " + "track secondaries:%s stack photons:%s", + m_maxNumPhotonsPerStep, m_maxBetaChangePerStep, + yes_no(m_trackSecondariesFirst), yes_no(m_stackPhotons)); + G4Cerenkov* process = new G4Cerenkov(name()); + process->SetVerboseLevel(m_verbosity); + process->SetMaxNumPhotonsPerStep(m_maxNumPhotonsPerStep); + process->SetMaxBetaChangePerStep(m_maxBetaChangePerStep); + process->SetTrackSecondariesFirst(m_trackSecondariesFirst); + process->SetStackPhotons(m_stackPhotons); + + auto pit = G4ParticleTable::GetParticleTable()->GetIterator(); + pit->reset(); + while( (*pit)() ){ + G4ParticleDefinition* particle = pit->value(); + if (process->IsApplicable(*particle)) { + G4ProcessManager* pmanager = particle->GetProcessManager(); + pmanager->AddProcess(process); + pmanager->SetProcessOrdering(process,idxPostStep); + } + } + } + + private: + double m_maxBetaChangePerStep; + int m_maxNumPhotonsPerStep; + int m_verbosity; + bool m_trackSecondariesFirst; + bool m_stackPhotons; + }; + } +} +#endif // DDG4_GEANT4CERENKOVPHYSICS_H + +#include "DDG4/Factories.h" +using namespace dd4hep::sim; +DECLARE_GEANT4ACTION(Geant4CerenkovPhysics) diff --git a/DDG4/plugins/Geant4DetectorGeometryConstruction.cpp b/DDG4/plugins/Geant4DetectorGeometryConstruction.cpp index 82eb1f5b4..51a2891d7 100644 --- a/DDG4/plugins/Geant4DetectorGeometryConstruction.cpp +++ b/DDG4/plugins/Geant4DetectorGeometryConstruction.cpp @@ -46,6 +46,8 @@ namespace dd4hep { bool m_debugPlacements = false; /// Property: Flag to debug regions during conversion mechanism bool m_debugRegions = false; + /// Property: Flag to debug regions during conversion mechanism + bool m_debugSurfaces = false; /// Property: Flag to dump all placements after the conversion procedure bool m_printPlacements = false; @@ -100,6 +102,7 @@ Geant4DetectorGeometryConstruction::Geant4DetectorGeometryConstruction(Geant4Con declareProperty("DebugVolumes", m_debugVolumes); declareProperty("DebugPlacements", m_debugPlacements); declareProperty("DebugRegions", m_debugRegions); + declareProperty("DebugSurfaces", m_debugSurfaces); declareProperty("PrintPlacements", m_printPlacements); declareProperty("PrintSensitives", m_printSensitives); @@ -126,6 +129,7 @@ void Geant4DetectorGeometryConstruction::constructGeo(Geant4DetectorConstruction conv.debugVolumes = m_debugVolumes; conv.debugPlacements = m_debugPlacements; conv.debugRegions = m_debugRegions; + conv.debugSurfaces = m_debugSurfaces; ctxt->geometry = conv.create(world).detach(); ctxt->geometry->printLevel = outputLevel(); diff --git a/DDG4/plugins/Geant4OpticalPhotonPhysics.cpp b/DDG4/plugins/Geant4OpticalPhotonPhysics.cpp new file mode 100644 index 000000000..76738aadf --- /dev/null +++ b/DDG4/plugins/Geant4OpticalPhotonPhysics.cpp @@ -0,0 +1,104 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// 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 +// +//========================================================================== + +/** \addtogroup Geant4PhysicsConstructor + * + * @{ + * \package Geant4Optical PhotonPhysics + * + * \brief PhysicsConstructor for Optical Photon physics + + This plugin allows to enable Optical Photon physics in the physics list + + * + * @} + */ +#ifndef DDG4_GEANT4OPTICALPHOTONPHYSICS_H +#define DDG4_GEANT4OPTICALPHOTONPHYSICS_H 1 + +// Framework include files +#include "DDG4/Geant4PhysicsList.h" + +/// Geant4 include files +#include "G4OpAbsorption.hh" +#include "G4OpRayleigh.hh" +#include "G4OpMieHG.hh" +#include "G4OpBoundaryProcess.hh" +#include "G4ParticleDefinition.hh" +#include "G4ParticleTypes.hh" +#include "G4ParticleTable.hh" +#include "G4ProcessManager.hh" + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace sim { + + /// Geant4 physics list action to enable OpticalPhoton physics + /** + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_SIMULATION + */ + class Geant4OpticalPhotonPhysics : public Geant4PhysicsList { + public: + /// Default constructor + Geant4OpticalPhotonPhysics() = delete; + /// Copy constructor + Geant4OpticalPhotonPhysics(const Geant4OpticalPhotonPhysics&) = delete; + /// Initializing constructor + Geant4OpticalPhotonPhysics(Geant4Context* ctxt, const std::string& nam) + : Geant4PhysicsList(ctxt, nam) + { + declareProperty("VerboseLevel", m_verbosity = 0); + } + /// Default destructor + virtual ~Geant4OpticalPhotonPhysics() = default; + /// Callback to construct processes (uses the G4 particle table) + virtual void constructProcesses(G4VUserPhysicsList* physics_list) { + this->Geant4PhysicsList::constructProcesses(physics_list); + info("+++ Constructing optical_photon processes:"); + info("+++ G4OpAbsorption G4OpRayleigh G4OpMieHG G4OpBoundaryProcess"); + G4ParticleTable* table = G4ParticleTable::GetParticleTable(); + G4ParticleDefinition* particle = table->FindParticle("opticalphoton"); + if (0 == particle) { + except("++ Cannot resolve 'opticalphoton' particle definition!"); + } + + G4OpBoundaryProcess* fBoundaryProcess = new G4OpBoundaryProcess(); + G4OpAbsorption* fAbsorptionProcess = new G4OpAbsorption(); + G4OpRayleigh* fRayleighScatteringProcess = new G4OpRayleigh(); + G4OpMieHG* fMieHGScatteringProcess = new G4OpMieHG(); + + fAbsorptionProcess->SetVerboseLevel(m_verbosity); + fRayleighScatteringProcess->SetVerboseLevel(m_verbosity); + fMieHGScatteringProcess->SetVerboseLevel(m_verbosity); + fBoundaryProcess->SetVerboseLevel(m_verbosity); + + G4ProcessManager* pmanager = particle->GetProcessManager(); + pmanager->AddDiscreteProcess(fAbsorptionProcess); + pmanager->AddDiscreteProcess(fRayleighScatteringProcess); + pmanager->AddDiscreteProcess(fMieHGScatteringProcess); + pmanager->AddDiscreteProcess(fBoundaryProcess); + } + private: + int m_verbosity; + }; + } +} +#endif // DDG4_GEANT4OPTICALPHOTONPHYSICS_H + +#include "DDG4/Factories.h" +using namespace dd4hep::sim; +DECLARE_GEANT4ACTION(Geant4OpticalPhotonPhysics) diff --git a/DDG4/plugins/Geant4ScintillationPhysics.cpp b/DDG4/plugins/Geant4ScintillationPhysics.cpp new file mode 100644 index 000000000..91b1eab0e --- /dev/null +++ b/DDG4/plugins/Geant4ScintillationPhysics.cpp @@ -0,0 +1,119 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// 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 +// +//========================================================================== + +/** \addtogroup Geant4PhysicsConstructor + * + * @{ + * \package Geant4ScintillationPhysics + * \brief PhysicsConstructor for Scintillation physics + + This plugin allows to enable Scintillation physics in the physics list + + * @} + */ +#ifndef DDG4_GEANT4SCINTILLATIONPHYSICS_H +#define DDG4_GEANT4SCINTILLATIONPHYSICS_H 1 + +/// Framework include files +#include "DDG4/Geant4PhysicsList.h" + +/// Geant4 include files +#include "G4ParticleTableIterator.hh" +#include "G4ParticleDefinition.hh" +#include "G4LossTableManager.hh" +#include "G4ProcessManager.hh" +#include "G4ParticleTypes.hh" +#include "G4ParticleTable.hh" +#include "G4EmSaturation.hh" +#include "G4Threading.hh" + +#include "G4Scintillation.hh" + +/// Namespace for the AIDA detector description toolkit +namespace dd4hep { + + /// Namespace for the Geant4 based simulation part of the AIDA detector description toolkit + namespace sim { + + /// Geant4 physics list action to enable Scintillation physics + /** + * \author M.Frank + * \version 1.0 + * \ingroup DD4HEP_SIMULATION + */ + class Geant4ScintillationPhysics : public Geant4PhysicsList { + public: + /// Default constructor + Geant4ScintillationPhysics() = delete; + /// Copy constructor + Geant4ScintillationPhysics(const Geant4ScintillationPhysics&) = delete; + /// Initializing constructor + Geant4ScintillationPhysics(Geant4Context* ctxt, const std::string& nam) + : Geant4PhysicsList(ctxt, nam) + { + declareProperty("ScintillationYieldFactor", m_scintillationYieldFactor = 1.0); + declareProperty("ScintillationExcitationRatio", m_scintillationExcitationRatio = 1.0); + declareProperty("FiniteRiseTime", m_finiteRiseTime = false); + declareProperty("TrackSecondariesFirst", m_trackSecondariesFirst = false); + declareProperty("StackPhotons", m_stackPhotons = true); + declareProperty("VerboseLevel", m_verbosity = 0); + } + /// Default destructor + virtual ~Geant4ScintillationPhysics() = default; + /// Callback to construct processes (uses the G4 particle table) + virtual void constructProcesses(G4VUserPhysicsList* physics_list) { + this->Geant4PhysicsList::constructProcesses(physics_list); + info("+++ Constructing: yield:%f excitation:%f rise-time:%s stack photons:%s", + m_scintillationYieldFactor, m_scintillationExcitationRatio, + yes_no(m_finiteRiseTime), yes_no(m_stackPhotons)); + G4Scintillation* process = new G4Scintillation(name()); + process->SetVerboseLevel(m_verbosity); + process->SetFiniteRiseTime(m_finiteRiseTime); + process->SetStackPhotons(m_stackPhotons); + process->SetTrackSecondariesFirst(m_trackSecondariesFirst); + process->SetScintillationYieldFactor(m_scintillationYieldFactor); + process->SetScintillationExcitationRatio(m_scintillationExcitationRatio); + // Use Birks Correction in the Scintillation process + if ( G4Threading::IsMasterThread() ) { + G4EmSaturation* emSaturation = + G4LossTableManager::Instance()->EmSaturation(); + process->AddSaturation(emSaturation); + } + + auto pit = G4ParticleTable::GetParticleTable()->GetIterator(); + pit->reset(); + while( (*pit)() ){ + G4ParticleDefinition* particle = pit->value(); + if (process->IsApplicable(*particle)) { + G4ProcessManager* pmanager = particle->GetProcessManager(); + pmanager->AddProcess(process); + pmanager->SetProcessOrdering(process,idxPostStep); + } + } + } + + private: + double m_scintillationYieldFactor; + double m_scintillationExcitationRatio; + int m_verbosity; + bool m_stackPhotons; + bool m_finiteRiseTime; + bool m_trackSecondariesFirst; + }; + } +} +#endif // DDG4_GEANT4SCINTILLATIONPHYSICS_H + +#include "DDG4/Factories.h" +using namespace dd4hep::sim; +DECLARE_GEANT4ACTION(Geant4ScintillationPhysics) diff --git a/DDG4/src/Geant4Converter.cpp b/DDG4/src/Geant4Converter.cpp index 3e4428f80..1b839ef23 100644 --- a/DDG4/src/Geant4Converter.cpp +++ b/DDG4/src/Geant4Converter.cpp @@ -401,31 +401,33 @@ void* Geant4Converter::handleMaterial(const string& name, Material medium) const mat = new G4Material(name, z, a, density, state, material->GetTemperature(), material->GetPressure()); } + stringstream str; + str << (*mat); + printout(lvl, "Geant4Converter", "++ Created G4 %s", str.str().c_str()); #if ROOT_VERSION_CODE > ROOT_VERSION(6,16,0) /// Attach the material properties if any G4MaterialPropertiesTable* tab = 0; TListIter propIt(&material->GetProperties()); - for(TObject* obj=*propIt; obj; obj = propIt.Next()) { - TNamed* n = (TNamed*)obj; - TGDMLMatrix *matrix = info.manager->GetGDMLMatrix(n->GetTitle()); - if ( 0 == tab ) { - tab = new G4MaterialPropertiesTable(); - mat->SetMaterialPropertiesTable(tab); - } + for(TObject* obj=propIt.Next(); obj; obj = propIt.Next()) { + TNamed* n = (TNamed*)obj; + TGDMLMatrix* matrix = info.manager->GetGDMLMatrix(n->GetTitle()); Geant4GeometryInfo::PropertyVector* v = (Geant4GeometryInfo::PropertyVector*)handleMaterialProperties(matrix); if ( 0 == v ) { except("Geant4Converter", "++ FAILED to create G4 material %s [Cannot convert property:%s]", material->GetName(), n->GetName()); } + if ( 0 == tab ) { + tab = new G4MaterialPropertiesTable(); + mat->SetMaterialPropertiesTable(tab); + } G4MaterialPropertyVector* vec = new G4MaterialPropertyVector(&v->bins[0], &v->values[0], v->bins.size()); + printout(lvl, "Geant4Converter", "++ Property: %-20s [%ld x %ld] -> %s ", + n->GetName(), matrix->GetRows(), matrix->GetCols(), n->GetTitle()); tab->AddProperty(n->GetName(), vec); } #endif - stringstream str; - str << (*mat); - printout(lvl, "Geant4Converter", "++ Created G4 %s", str.str().c_str()); } info.g4Materials[medium] = mat; } @@ -1060,6 +1062,7 @@ void* Geant4Converter::handleMaterialProperties(TObject* matrix) const { Geant4GeometryInfo& info = data(); Geant4GeometryInfo::PropertyVector* g4 = info.g4OpticalProperties[m]; if (!g4) { + PrintLevel lvl = debugMaterials ? ALWAYS : outputLevel; g4 = new Geant4GeometryInfo::PropertyVector(); size_t rows = m->GetRows(); g4->name = m->GetName(); @@ -1067,9 +1070,11 @@ void* Geant4Converter::handleMaterialProperties(TObject* matrix) const { g4->bins.reserve(rows); g4->values.reserve(rows); for(size_t i=0; i<rows; ++i) { - g4->bins.push_back(m->Get(0,i) /* *CLHEP::eV/units::eV */); - g4->values.push_back(m->Get(1,i)); + g4->bins.push_back(m->Get(i,0) /* *CLHEP::eV/units::eV */); + g4->values.push_back(m->Get(i,1)); } + printout(lvl, "Geant4Converter", "++ Successfully converted material property:%s : %s [%ld rows]", + m->GetName(), m->GetTitle(), rows); info.g4OpticalProperties[m] = g4; } return g4; @@ -1177,9 +1182,15 @@ void* Geant4Converter::handleOpticalSurface(TObject* surface) const { g4 = new G4OpticalSurface(s->GetName(), model, finish, type, s->GetValue()); g4->SetSigmaAlpha(s->GetSigmaAlpha()); // not implemented: g4->SetPolish(s->GetPolish()); + printout(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter", + "++ Created OpticalSurface: %-18s type:%s model:%s finish:%s", + s->GetName(), + TGeoOpticalSurface::TypeToString(s->GetType()), + TGeoOpticalSurface::ModelToString(s->GetModel()), + TGeoOpticalSurface::FinishToString(s->GetFinish())); G4MaterialPropertiesTable* tab = 0; TListIter it(&s->GetProperties()); - for(TObject* obj = *it; obj; obj = it.Next()) { + for(TObject* obj = it.Next(); obj; obj = it.Next()) { TNamed* n = (TNamed*)obj; TGDMLMatrix *matrix = info.manager->GetGDMLMatrix(n->GetTitle()); if ( 0 == tab ) { @@ -1195,6 +1206,9 @@ void* Geant4Converter::handleOpticalSurface(TObject* surface) const { G4MaterialPropertyVector* vec = new G4MaterialPropertyVector(&v->bins[0], &v->values[0], v->bins.size()); tab->AddProperty(n->GetName(), vec); + printout(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter", + "++ Property: %-20s [%ld x %ld] --> %s", + n->GetName(), matrix->GetRows(), matrix->GetCols(), n->GetTitle()); } info.g4OpticalSurfaces[s] = g4; } @@ -1210,6 +1224,9 @@ void* Geant4Converter::handleSkinSurface(TObject* surface) const { G4OpticalSurface* s = info.g4OpticalSurfaces[OpticalSurface(surf->GetSurface())]; G4LogicalVolume* v = info.g4Volumes[surf->GetVolume()]; g4 = new G4LogicalSkinSurface(surf->GetName(), v, s); + printout(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter", + "++ Created SkinSurface: %-18s optical:%s", + surf->GetName(), surf->GetSurface()->GetName()); info.g4SkinSurfaces[surf] = g4; } return g4; @@ -1225,6 +1242,9 @@ void* Geant4Converter::handleBorderSurface(TObject* surface) const { G4VPhysicalVolume* n1 = info.g4Placements[surf->GetNode1()]; G4VPhysicalVolume* n2 = info.g4Placements[surf->GetNode2()]; g4 = new G4LogicalBorderSurface(surf->GetName(), n1, n2, s); + printout(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter", + "++ Created BorderSurface: %-18s optical:%s", + surf->GetName(), surf->GetSurface()->GetName()); info.g4BorderSurfaces[surf] = g4; } return g4; diff --git a/DDG4/src/Geant4PhysicsList.cpp b/DDG4/src/Geant4PhysicsList.cpp index c07fa6a3d..2fd52f920 100644 --- a/DDG4/src/Geant4PhysicsList.cpp +++ b/DDG4/src/Geant4PhysicsList.cpp @@ -44,17 +44,15 @@ namespace { }; struct EmptyPhysics : public G4VModularPhysicsList { - EmptyPhysics() {} - virtual ~EmptyPhysics() {} - virtual void ConstructProcess() {} - virtual void ConstructParticle() {} + EmptyPhysics() {} + virtual ~EmptyPhysics() {} }; struct ParticlePhysics : public G4VPhysicsConstructor { Geant4PhysicsListActionSequence* seq; G4VUserPhysicsList* phys; ParticlePhysics(Geant4PhysicsListActionSequence* s, G4VUserPhysicsList* p) - : seq(s), phys(p) {} - virtual ~ParticlePhysics() {} + : seq(s), phys(p) {} + virtual ~ParticlePhysics() {} virtual void ConstructProcess() { seq->constructProcesses(phys); if ( seq->transportation() ) { @@ -62,21 +60,23 @@ namespace { ph->AddTransportation(); } } - virtual void ConstructParticle() { seq->constructParticles(phys); } + virtual void ConstructParticle() { + seq->constructParticles(phys); + } }; } /// Default constructor Geant4PhysicsList::Process::Process() - : ordAtRestDoIt(-1), ordAlongSteptDoIt(-1), ordPostStepDoIt(-1) { + : ordAtRestDoIt(-1), ordAlongSteptDoIt(-1), ordPostStepDoIt(-1) { } /// Copy constructor Geant4PhysicsList::Process::Process(const Process& p) - : name(p.name), ordAtRestDoIt(p.ordAtRestDoIt), ordAlongSteptDoIt(p.ordAlongSteptDoIt), ordPostStepDoIt(p.ordPostStepDoIt) { + : name(p.name), ordAtRestDoIt(p.ordAtRestDoIt), ordAlongSteptDoIt(p.ordAlongSteptDoIt), ordPostStepDoIt(p.ordPostStepDoIt) { } /// Assignment operator -Geant4PhysicsList::Process& Geant4PhysicsList::Process::operator=(const Process& p) { +Geant4PhysicsList::Process& Geant4PhysicsList::Process::operator=(const Process& p) { if ( this != &p ) { name = p.name; ordAtRestDoIt = p.ordAtRestDoIt; @@ -88,12 +88,12 @@ Geant4PhysicsList::Process& Geant4PhysicsList::Process::operator=(const Process& /// Standard constructor Geant4PhysicsList::Geant4PhysicsList(Geant4Context* ctxt, const string& nam) - : Geant4Action(ctxt, nam) { + : Geant4Action(ctxt, nam) { InstanceCount::increment(this); } /// Default destructor -Geant4PhysicsList::~Geant4PhysicsList() { +Geant4PhysicsList::~Geant4PhysicsList() { InstanceCount::decrement(this); } @@ -111,11 +111,21 @@ void Geant4PhysicsList::dump() { printout(ALWAYS,name(),"+++ PhysicsConstructor: %s",(*i).c_str()); for (ParticleConstructors::const_iterator i = m_particles.begin(); i != m_particles.end(); ++i) printout(ALWAYS,name(),"+++ ParticleConstructor: %s",(*i).c_str()); - for (PhysicsProcesses::const_iterator i = m_processes.begin(); i != m_processes.end(); ++i) { + for (ParticleConstructors::const_iterator i = m_particlegroups.begin(); i != m_particlegroups.end(); ++i) + printout(ALWAYS,name(),"+++ ParticleGroupConstructor: %s",(*i).c_str()); + for (PhysicsProcesses::const_iterator i = m_discreteProcesses.begin(); i != m_discreteProcesses.end(); ++i) { + const string& part_name = (*i).first; + const ParticleProcesses& procs = (*i).second; + printout(ALWAYS,name(),"+++ DiscretePhysicsProcesses of particle %s",part_name.c_str()); + for (ParticleProcesses::const_iterator ip = procs.begin(); ip != procs.end(); ++ip) { + printout(ALWAYS,name(),"+++ Process %s", (*ip).name.c_str()); + } + } + for (PhysicsProcesses::const_iterator i = m_processes.begin(); i != m_processes.end(); ++i) { const string& part_name = (*i).first; const ParticleProcesses& procs = (*i).second; printout(ALWAYS,name(),"+++ PhysicsProcesses of particle %s",part_name.c_str()); - for (ParticleProcesses::const_iterator ip = procs.begin(); ip != procs.end(); ++ip) { + for (ParticleProcesses::const_iterator ip = procs.begin(); ip != procs.end(); ++ip) { const Process& p = (*ip); printout(ALWAYS,name(),"+++ Process %s ordAtRestDoIt=%d ordAlongSteptDoIt=%d ordPostStepDoIt=%d", p.name.c_str(),p.ordAtRestDoIt,p.ordAlongSteptDoIt,p.ordPostStepDoIt); @@ -148,15 +158,24 @@ void Geant4PhysicsList::addParticleProcess(const std::string& part_name, processes(part_name).push_back(p); } +/// Add discrete particle process by name with arguments +void Geant4PhysicsList::addDiscreteParticleProcess(const std::string& part_name, + const std::string& proc_name) +{ + Process p; + p.name = proc_name; + discreteProcesses(part_name).push_back(p); +} + /// Add PhysicsConstructor by name void Geant4PhysicsList::addPhysicsConstructor(const std::string& phys_name) { physics().push_back(phys_name); } /// Access processes for one particle type -Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::processes(const string& nam) { +Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::processes(const string& nam) { PhysicsProcesses::iterator i = m_processes.find(nam); - if (i != m_processes.end()) { + if (i != m_processes.end()) { return (*i).second; } pair<PhysicsProcesses::iterator, bool> ret = m_processes.insert(make_pair(nam, ParticleProcesses())); @@ -166,7 +185,27 @@ Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::processes(const string& /// Access processes for one particle type (CONST) const Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::processes(const string& nam) const { PhysicsProcesses::const_iterator i = m_processes.find(nam); - if (i != m_processes.end()) { + if (i != m_processes.end()) { + return (*i).second; + } + except("Failed to access the physics process '%s' [Unknown-Process]", nam.c_str()); + throw runtime_error("Failed to access the physics process"); // never called anyway +} + +/// Access discrete processes for one particle type +Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::discreteProcesses(const string& nam) { + PhysicsProcesses::iterator i = m_discreteProcesses.find(nam); + if (i != m_discreteProcesses.end()) { + return (*i).second; + } + pair<PhysicsProcesses::iterator, bool> ret = m_discreteProcesses.insert(make_pair(nam, ParticleProcesses())); + return (*(ret.first)).second; +} + +/// Access discrete processes for one particle type (CONST) +const Geant4PhysicsList::ParticleProcesses& Geant4PhysicsList::discreteProcesses(const string& nam) const { + PhysicsProcesses::const_iterator i = m_discreteProcesses.find(nam); + if (i != m_discreteProcesses.end()) { return (*i).second; } except("Failed to access the physics process '%s' [Unknown-Process]", nam.c_str()); @@ -190,13 +229,13 @@ void Geant4PhysicsList::adoptPhysicsConstructor(Geant4Action* action) { } /// Callback to construct particle decays -void Geant4PhysicsList::constructPhysics(G4VModularPhysicsList* physics_pointer) { +void Geant4PhysicsList::constructPhysics(G4VModularPhysicsList* physics_pointer) { debug("constructPhysics %p", physics_pointer); - for (PhysicsConstructors::iterator i=m_physics.begin(); i != m_physics.end(); ++i) { + for (PhysicsConstructors::iterator i=m_physics.begin(); i != m_physics.end(); ++i) { PhysicsConstructors::value_type& ctor = *i; if ( 0 == ctor.pointer ) { G4VPhysicsConstructor* p = PluginService::Create<G4VPhysicsConstructor*>(ctor); - if (!p) { + if (!p) { except("Failed to create the physics entities " "for the G4VPhysicsConstructor '%s'", ctor.c_str()); } @@ -208,49 +247,73 @@ void Geant4PhysicsList::constructPhysics(G4VModularPhysicsList* physics_pointer) } /// constructParticle callback -void Geant4PhysicsList::constructParticles(G4VUserPhysicsList* physics_pointer) { +void Geant4PhysicsList::constructParticles(G4VUserPhysicsList* physics_pointer) { debug("constructParticles %p", physics_pointer); /// Now define all particles - for (ParticleConstructors::const_iterator i = m_particles.begin(); i != m_particles.end(); ++i) { + for (ParticleConstructors::const_iterator i = m_particles.begin(); i != m_particles.end(); ++i) { const ParticleConstructors::value_type& ctor = *i; G4ParticleDefinition* def = PluginService::Create<G4ParticleDefinition*>(ctor); - if (!def) { + if (!def) { /// Check if we have here a particle group constructor long* result = (long*) PluginService::Create<long>(ctor); - if (!result || *result != 1L) { + if (!result || *result != 1L) { except("Failed to create particle type '%s' result=%d", ctor.c_str(), result); } info("Constructed Geant4 particle %s [using signature long (*)()]",ctor.c_str()); } } /// Now define all particles - for (ParticleConstructors::const_iterator i = m_particlegroups.begin(); i != m_particlegroups.end(); ++i) { + for (ParticleConstructors::const_iterator i = m_particlegroups.begin(); i != m_particlegroups.end(); ++i) { const ParticleConstructors::value_type& ctor = *i; /// Check if we have here a particle group constructor long* result = (long*) PluginService::Create<long>(ctor); - if (!result || *result != 1L) { + if (!result || *result != 1L) { except("Failed to create particle type '%s' result=%d", ctor.c_str(), result); } + info("Constructed Geant4 particle group %s [using signature long (*)()]",ctor.c_str()); } } /// Callback to construct processes (uses the G4 particle table) -void Geant4PhysicsList::constructProcesses(G4VUserPhysicsList* physics_pointer) { +void Geant4PhysicsList::constructProcesses(G4VUserPhysicsList* physics_pointer) { debug("constructProcesses %p", physics_pointer); - for (PhysicsProcesses::const_iterator i = m_processes.begin(); i != m_processes.end(); ++i) { + for (PhysicsProcesses::const_iterator i = m_discreteProcesses.begin(); i != m_discreteProcesses.end(); ++i) { + const string& part_name = (*i).first; + const ParticleProcesses& procs = (*i).second; + vector<G4ParticleDefinition*> defs(Geant4ParticleHandle::g4DefinitionsRegEx(part_name)); + if (defs.empty()) { + except("Particle:%s Cannot find the corresponding entry in the particle table.", part_name.c_str()); + } + for (vector<G4ParticleDefinition*>::const_iterator id = defs.begin(); id != defs.end(); ++id) { + G4ParticleDefinition* particle = *id; + G4ProcessManager* mgr = particle->GetProcessManager(); + for (ParticleProcesses::const_iterator ip = procs.begin(); ip != procs.end(); ++ip) { + const Process& p = (*ip); + G4VProcess* g4 = PluginService::Create<G4VProcess*>(p.name); + if (!g4) { // Error no factory for this process + except("Particle:%s -> [%s] Cannot create discrete physics process %s", + part_name.c_str(), particle->GetParticleName().c_str(), p.name.c_str()); + } + mgr->AddDiscreteProcess(g4); + info("Particle:%s -> [%s] added discrete process %s", + part_name.c_str(), particle->GetParticleName().c_str(), p.name.c_str()); + } + } + } + for (PhysicsProcesses::const_iterator i = m_processes.begin(); i != m_processes.end(); ++i) { const string& part_name = (*i).first; const ParticleProcesses& procs = (*i).second; vector<G4ParticleDefinition*> defs(Geant4ParticleHandle::g4DefinitionsRegEx(part_name)); - if (defs.empty()) { + if (defs.empty()) { except("Particle:%s Cannot find the corresponding entry in the particle table.", part_name.c_str()); } - for (vector<G4ParticleDefinition*>::const_iterator id = defs.begin(); id != defs.end(); ++id) { + for (vector<G4ParticleDefinition*>::const_iterator id = defs.begin(); id != defs.end(); ++id) { G4ParticleDefinition* particle = *id; G4ProcessManager* mgr = particle->GetProcessManager(); - for (ParticleProcesses::const_iterator ip = procs.begin(); ip != procs.end(); ++ip) { + for (ParticleProcesses::const_iterator ip = procs.begin(); ip != procs.end(); ++ip) { const Process& p = (*ip); G4VProcess* g4 = PluginService::Create<G4VProcess*>(p.name); - if (!g4) { // Error no factory for this process + if (!g4) { // Error no factory for this process except("Particle:%s -> [%s] Cannot create physics process %s", part_name.c_str(), particle->GetParticleName().c_str(), p.name.c_str()); } @@ -269,7 +332,7 @@ void Geant4PhysicsList::enable(G4VUserPhysicsList* /* physics */) { /// Standard constructor Geant4PhysicsListActionSequence::Geant4PhysicsListActionSequence(Geant4Context* ctxt, const string& nam) - : Geant4Action(ctxt, nam), m_transportation(false), m_decays(false), m_rangecut(0.7*CLHEP::mm) { + : Geant4Action(ctxt, nam), m_transportation(false), m_decays(false), m_rangecut(0.7*CLHEP::mm) { declareProperty("transportation", m_transportation); declareProperty("extends", m_extends); declareProperty("decays", m_decays); @@ -279,7 +342,7 @@ Geant4PhysicsListActionSequence::Geant4PhysicsListActionSequence(Geant4Context* } /// Default destructor -Geant4PhysicsListActionSequence::~Geant4PhysicsListActionSequence() { +Geant4PhysicsListActionSequence::~Geant4PhysicsListActionSequence() { m_actors(&Geant4Action::release); m_actors.clear(); m_process.clear(); @@ -308,7 +371,7 @@ void Geant4PhysicsListActionSequence::installCommandMessenger() { /// Dump content to stdout void Geant4PhysicsListActionSequence::dump() { - printout(ALWAYS,name(),"+++ Dump"); + printout(ALWAYS,name(),"+++ Dump of physics list component(s)"); printout(ALWAYS,name(),"+++ Extension name %s",m_extends.c_str()); printout(ALWAYS,name(),"+++ Transportation flag: %d",m_transportation); printout(ALWAYS,name(),"+++ Program decays: %d",m_decays); @@ -317,8 +380,8 @@ void Geant4PhysicsListActionSequence::dump() { } /// Add an actor responding to all callbacks. Sequence takes ownership. -void Geant4PhysicsListActionSequence::adopt(Geant4PhysicsList* action) { - if (action) { +void Geant4PhysicsListActionSequence::adopt(Geant4PhysicsList* action) { + if (action) { action->addRef(); m_actors.add(action); return; @@ -327,38 +390,38 @@ void Geant4PhysicsListActionSequence::adopt(Geant4PhysicsList* action) { } /// Callback to construct particles -void Geant4PhysicsListActionSequence::constructParticles(G4VUserPhysicsList* physics_pointer) { +void Geant4PhysicsListActionSequence::constructParticles(G4VUserPhysicsList* physics_pointer) { m_particle(physics_pointer); m_actors(&Geant4PhysicsList::constructParticles, physics_pointer); } /// Callback to execute physics constructors -void Geant4PhysicsListActionSequence::constructPhysics(G4VModularPhysicsList* physics_pointer) { +void Geant4PhysicsListActionSequence::constructPhysics(G4VModularPhysicsList* physics_pointer) { m_physics(physics_pointer); m_actors(&Geant4PhysicsList::constructPhysics, physics_pointer); } /// constructProcess callback -void Geant4PhysicsListActionSequence::constructProcesses(G4VUserPhysicsList* physics_pointer) { +void Geant4PhysicsListActionSequence::constructProcesses(G4VUserPhysicsList* physics_pointer) { m_actors(&Geant4PhysicsList::constructProcesses, physics_pointer); m_process(physics_pointer); - if (m_decays) { + if (m_decays) { constructDecays(physics_pointer); } } /// Callback to construct particle decays -void Geant4PhysicsListActionSequence::constructDecays(G4VUserPhysicsList* physics_pointer) { +void Geant4PhysicsListActionSequence::constructDecays(G4VUserPhysicsList* physics_pointer) { G4ParticleTable* pt = G4ParticleTable::GetParticleTable(); G4ParticleTable::G4PTblDicIterator* iter = pt->GetIterator(); // Add Decay Process G4Decay* decay = new G4Decay(); info("ConstructDecays %p",physics_pointer); iter->reset(); - while ((*iter)()) { + while ((*iter)()) { G4ParticleDefinition* p = iter->value(); G4ProcessManager* mgr = p->GetProcessManager(); - if (decay->IsApplicable(*p)) { + if (decay->IsApplicable(*p)) { mgr->AddProcess(decay); // set ordering for PostStepDoIt and AtRestDoIt mgr->SetProcessOrdering(decay, idxPostStep); diff --git a/examples/OpticalSurfaces/CMakeLists.txt b/examples/OpticalSurfaces/CMakeLists.txt index 3a57a77de..1f51c6910 100644 --- a/examples/OpticalSurfaces/CMakeLists.txt +++ b/examples/OpticalSurfaces/CMakeLists.txt @@ -20,15 +20,15 @@ dd4hep_package ( OpticalSurfaces MAJOR 0 MINOR 0 PATCH 1 INCLUDE_DIRS include ) set(OpticalSurfaces_INSTALL ${CMAKE_INSTALL_PREFIX}/examples/OpticalSurfaces) -# dd4hep_add_plugin( OpticalSurfacesExample SOURCES src/*.cpp ) -dd4hep_install_dir( compact DESTINATION ${OpticalSurfaces_INSTALL} ) +dd4hep_add_plugin( OpticalSurfacesExample SOURCES src/*.cpp ) +dd4hep_install_dir( compact scripts DESTINATION ${OpticalSurfaces_INSTALL} ) dd4hep_configure_scripts( OpticalSurfaces DEFAULT_SETUP WITH_TESTS) # message (STATUS "OpticalSurfaces: ROOT version: ${ROOT_VERSION}") # #---Testing: Load ROOT GDMLMatrix objects from compact # and register them to the TGeoManager -dd4hep_add_test_reg( Surfaces_read_GDMLMatrix +dd4hep_add_test_reg( Surfaces_read_GDMLMatrices COMMAND "${CMAKE_INSTALL_PREFIX}/bin/run_test_OpticalSurfaces.sh" EXEC_ARGS geoPluginRun -volmgr -destroy -compact file:${OpticalSurfaces_INSTALL}/compact/ReadGDMLMatrices.xml @@ -49,3 +49,13 @@ dd4hep_add_test_reg( Surfaces_read_OpticalSurfaces REGEX_FAIL " ERROR ;EXCEPTION;Exception" ) # +#---Testing: Load ROOT material properties and dump material +dd4hep_add_test_reg( Surfaces_read_MaterialProperties + COMMAND "${CMAKE_INSTALL_PREFIX}/bin/run_test_OpticalSurfaces.sh" + EXEC_ARGS geoPluginRun -volmgr -destroy -print INFO + -compact file:${OpticalSurfaces_INSTALL}/compact/ReadMaterialProperties.xml + -plugin DD4hep_MaterialTable -name Water + REGEX_PASS "Property: SLOWCOMPONENT \\[32 x 2\\] --> SLOWCOMPONENT__0x123aff00" + REGEX_FAIL " ERROR ;EXCEPTION;Exception" + ) +# diff --git a/examples/OpticalSurfaces/compact/OpNovice.xml b/examples/OpticalSurfaces/compact/OpNovice.xml new file mode 100644 index 000000000..d9438a697 --- /dev/null +++ b/examples/OpticalSurfaces/compact/OpNovice.xml @@ -0,0 +1,375 @@ +<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="TestMaterialProperties" + title="Test reading of TGeo's Material Properties" + author="Markus Frank" + url="None" + status="development" + version="1.0"> + <comment>Test reading of TGeo's Material Properties</comment> + </info> + + <define> + <constant name="world_side" value="1*m"/> + <constant name="world_x" value="world_side/2"/> + <constant name="world_y" value="world_side/2"/> + <constant name="world_z" value="world_side/2"/> + </define> + <debug> + <type name="surface" value="1"/> + </debug> + + <properties> + <matrix name="RINDEX__Air" coldim="2" values=" + 2.034*eV 1. + 2.068*eV 1. + 2.103*eV 1. + 2.139*eV 1. + 2.177*eV 1. + 2.216*eV 1. + 2.256*eV 1. + 2.298*eV 1. + 2.341*eV 1. + 2.386*eV 1. + 2.433*eV 1. + 2.481*eV 1. + 2.532*eV 1. + 2.585*eV 1. + 2.640*eV 1. + 2.697*eV 1. + 2.757*eV 1. + 2.820*eV 1. + 2.885*eV 1. + 2.954*eV 1. + 3.026*eV 1. + 3.102*eV 1. + 3.181*eV 1. + 3.265*eV 1. + 3.353*eV 1. + 3.446*eV 1. + 3.545*eV 1. + 3.649*eV 1. + 3.760*eV 1. + 3.877*eV 1. + 4.002*eV 1. + 4.136*eV 1. + "/> + <matrix name="RINDEX__Water" coldim="2" values=" + 2.034*eV 1.3435 + 2.068*eV 1.344 + 2.103*eV 1.3445 + 2.139*eV 1.345 + 2.177*eV 1.3455 + 2.216*eV 1.346 + 2.256*eV 1.3465 + 2.298*eV 1.347 + 2.341*eV 1.3475 + 2.386*eV 1.348 + 2.433*eV 1.3485 + 2.481*eV 1.3492 + 2.532*eV 1.35 + 2.585*eV 1.3505 + 2.640*eV 1.351 + 2.697*eV 1.3518 + 2.757*eV 1.3522 + 2.820*eV 1.3530 + 2.885*eV 1.3535 + 2.954*eV 1.354 + 3.026*eV 1.3545 + 3.102*eV 1.355 + 3.181*eV 1.3555 + 3.265*eV 1.356 + 3.353*eV 1.3568 + 3.446*eV 1.3572 + 3.545*eV 1.358 + 3.649*eV 1.3585 + 3.760*eV 1.359 + 3.877*eV 1.3595 + 4.002*eV 1.36 + 4.136*eV 1.3608 + "/> + <matrix name="ABSLENGTH__Water" coldim="2" values=" + 2.034*eV 3.448*m + 2.068*eV 4.082*m + 2.103*eV 6.329*m + 2.139*eV 9.174*m + 2.177*eV 12.346*m + 2.216*eV 13.889*m + 2.256*eV 15.152*m + 2.298*eV 17.241*m + 2.341*eV 18.868*m + 2.386*eV 20.000*m + 2.433*eV 26.316*m + 2.481*eV 35.714*m + 2.532*eV 45.455*m + 2.585*eV 47.619*m + 2.640*eV 52.632*m + 2.697*eV 52.632*m + 2.757*eV 55.556*m + 2.820*eV 52.632*m + 2.885*eV 52.632*m + 2.954*eV 47.619*m + 3.026*eV 45.455*m + 3.102*eV 41.667*m + 3.181*eV 37.037*m + 3.265*eV 33.333*m + 3.353*eV 30.000*m + 3.446*eV 28.500*m + 3.545*eV 27.000*m + 3.649*eV 24.500*m + 3.760*eV 22.000*m + 3.877*eV 19.500*m + 4.002*eV 17.500*m + 4.136*eV 14.500*m + "/> + <matrix name= "FASTCOMPONENT__Water" coldim="2" values=" + 2.034*eV 1 + 2.068*eV 1 + 2.103*eV 1 + 2.139*eV 1 + 2.177*eV 1 + 2.216*eV 1 + 2.256*eV 1 + 2.298*eV 1 + 2.341*eV 1 + 2.386*eV 1 + 2.433*eV 1 + 2.481*eV 1 + 2.532*eV 1 + 2.585*eV 1 + 2.640*eV 1 + 2.697*eV 1 + 2.757*eV 1 + 2.820*eV 1 + 2.885*eV 1 + 2.954*eV 1 + 3.026*eV 1 + 3.102*eV 1 + 3.181*eV 1 + 3.265*eV 1 + 3.353*eV 1 + 3.446*eV 1 + 3.545*eV 1 + 3.649*eV 1 + 3.760*eV 1 + 3.877*eV 1 + 4.002*eV 1 + 4.136*eV 1 + "/> + <matrix name= "SLOWCOMPONENT__Water" coldim="2" values=" + 2.034*eV 0.01 + 2.068*eV 1 + 2.103*eV 2 + 2.139*eV 3 + 2.177*eV 4 + 2.216*eV 5 + 2.256*eV 6 + 2.298*eV 7 + 2.341*eV 8 + 2.386*eV 9 + 2.433*eV 8 + 2.481*eV 7 + 2.532*eV 6 + 2.585*eV 4 + 2.640*eV 3 + 2.697*eV 2 + 2.757*eV 1 + 2.820*eV 0.01 + 2.885*eV 1 + 2.954*eV 2 + 3.026*eV 3 + 3.102*eV 4 + 3.181*eV 5 + 3.265*eV 6 + 3.353*eV 7 + 3.446*eV 8 + 3.545*eV 9 + 3.649*eV 8 + 3.760*eV 7 + 3.877*eV 6 + 4.002*eV 5 + 4.136*eV 4 + "/> + <matrix name= "MIEHG__Water" coldim="2" values=" + 1.56962*eV 167024.4*m + 1.58974*eV 158726.7*m + 1.61039*eV 150742.0*m + 1.63157*eV 143062.5*m + 1.65333*eV 135680.2*m + 1.67567*eV 128587.4*m + 1.69863*eV 121776.3*m + 1.72222*eV 115239.5*m + 1.74647*eV 108969.5*m + 1.77142*eV 102958.8*m + 1.79710*eV 97200.35*m + 1.82352*eV 91686.86*m + 1.85074*eV 86411.33*m + 1.87878*eV 81366.79*m + 1.90769*eV 76546.42*m + 1.93749*eV 71943.46*m + 1.96825*eV 67551.29*m + 1.99999*eV 63363.36*m + 2.03278*eV 59373.25*m + 2.06666*eV 55574.61*m + 2.10169*eV 51961.24*m + 2.13793*eV 48527.00*m + 2.17543*eV 45265.87*m + 2.21428*eV 42171.94*m + 2.25454*eV 39239.39*m + 2.29629*eV 36462.50*m + 2.33962*eV 33835.68*m + 2.38461*eV 31353.41*m + 2.43137*eV 29010.30*m + 2.47999*eV 26801.03*m + 2.53061*eV 24720.42*m + 2.58333*eV 22763.36*m + 2.63829*eV 20924.88*m + 2.69565*eV 19200.07*m + 2.75555*eV 17584.16*m + 2.81817*eV 16072.45*m + 2.88371*eV 14660.38*m + 2.95237*eV 13343.46*m + 3.02438*eV 12117.33*m + 3.09999*eV 10977.70*m + 3.17948*eV 9920.416*m + 3.26315*eV 8941.407*m + 3.35134*eV 8036.711*m + 3.44444*eV 7202.470*m + 3.54285*eV 6434.927*m + 3.64705*eV 5730.429*m + 3.75757*eV 5085.425*m + 3.87499*eV 4496.467*m + 3.99999*eV 3960.210*m + 4.13332*eV 3473.413*m + 4.27585*eV 3032.937*m + 4.42856*eV 2635.746*m + 4.59258*eV 2278.907*m + 4.76922*eV 1959.588*m + 4.95999*eV 1675.064*m + 5.16665*eV 1422.710*m + 5.39129*eV 1200.004*m + 5.63635*eV 1004.528*m + 5.90475*eV 833.9666*m + 6.19998*eV 686.1063*m + "/> + <ignore> + <matrix name= "REFLECTIVITY__Water" coldim="2" values=" + 2.034*eV + 2.068*eV + 2.103*eV + 2.139*eV + 2.177*eV + 2.216*eV + 2.256*eV + 2.298*eV + 2.341*eV + 2.386*eV + 2.433*eV + 2.481*eV + 2.532*eV + 2.585*eV + 2.640*eV + 2.697*eV + 2.757*eV + 2.820*eV + 2.885*eV + 2.954*eV + 3.026*eV + 3.102*eV + 3.181*eV + 3.265*eV + 3.353*eV + 3.446*eV + 3.545*eV + 3.649*eV + 3.760*eV + 3.877*eV + 4.002*eV + 4.136*eV + "/> + </ignore> + </properties> + + <materials> + <element Z="1" formula="H" name="H" > + <atom type="A" unit="g/mol" value="1.00794" /> + </element> + <element Z="8" formula="O" name="O" > + <atom type="A" unit="g/mol" value="15.9994" /> + </element> + <element Z="7" formula="N" name="N" > + <atom type="A" unit="g/mol" value="14.0068" /> + </element> + <element Z="18" formula="Ar" name="Ar" > + <atom type="A" unit="g/mol" value="39.9477" /> + </element> + + <material name="Air"> + <D type="density" unit="g/cm3" value="0.0012"/> + <fraction n="0.754" ref="N"/> + <fraction n="0.234" ref="O"/> + <fraction n="0.012" ref="Ar"/> + <property name="RINDEX" ref="RINDEX__Air"/> + </material> + <!-- We model vakuum just as very thin air --> + <material name="Vacuum"> + <D type="density" unit="g/cm3" value="0.0000000001" /> + <fraction n="0.754" ref="N"/> + <fraction n="0.234" ref="O"/> + <fraction n="0.012" ref="Ar"/> + </material> + <material name="Water"> + <D type="density" value="1.0" unit="g/cm3"/> + <composite n="2" ref="H"/> + <composite n="1" ref="O"/> + <property name="RINDEX" ref="RINDEX__Water"/> + <property name="ABSLENGTH" ref="ABSLENGTH__Water"/> + <property name="FASTCOMPONENT" ref="FASTCOMPONENT__Water"/> + <property name="SLOWCOMPONENT" ref="SLOWCOMPONENT__Water"/> + <property name="MIEHG" ref="MIEHG__Water"/> + <const_property name="SCINTILLATIONYIELD" value="50.0/MeV"/> + <const_property name="RESOLUTIONSCALE" value="1.0"/> + <const_property name="FASTTIMECONSTANT" value="1.0*ns"/> + <const_property name="SLOWTIMECONSTANT" value="10.0*ns"/> + <const_property name="YIELDRATIO" value="0.8"/> + <const_property name="MIEHG_FORWARD" value="0.99"/> + <const_property name="MIEHG_BACKWARD" value="0.99"/> + <const_property name="MIEHG_FORWARD_RATIO" value="0.8"/> + </material> + </materials> + + <display> + <vis name="TankVis" alpha="0.5" r="0" g="0.0" b="1.0" showDaughters="true" visible="true"/> + <vis name="BubbleVis" alpha="0.3" r="1" g="0.0" b="0.0" showDaughters="true" visible="true"/> + </display> + + <surfaces> + <alt___opticalsurface name="/world/BubbleDevice#WaterSurface" finish="ground" model="unified" type="dielectric_dielectric"> + <property name="RINDEX" coldim="2" values="2.034*eV 1.35 4.136*eV 1.40"/> + <property name="SPECULARLOBECONSTANT" coldim="2" values="2.034*eV 0.3 4.136*eV 0.3 "/> + <property name="SPECULARSPIKECONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/> + <property name="BACKSCATTERCONSTANT" coldim="2" values="2.034*eV 0.2 4.136*eV 0.2 "/> + </alt___opticalsurface> + <opticalsurface name="/world/BubbleDevice#WaterSurface" finish="Rough_LUT" model="DAVIS" type="dielectric_LUTDAVIS"/> + <opticalsurface name="/world/BubbleDevice#AirSurface" finish="polished" model="glisur" type="dielectric_dielectric"> + <property name="REFLECTIVITY" coldim="2" values="2.034*eV 0.3 4.136*eV 0.5"/> + <property name="EFFICIENCY" coldim="2" values="2.034*eV 0.8 4.136*eV 1.0"/> + </opticalsurface> + </surfaces> + + <detectors> + <detector name="BubbleDevice" type="DD4hep_OpNovice" vis="TankVis" id ="1" sensitive="true" readout="BubbleDeviceHits"> + <tank x="1*m" y="1*m" z="1*m" material="Water" vis="TankVis"/> + <bubble x="0.5*m" y="0.5*m" z="0.5*m" material="Air" vis="BubbleVis"/> + </detector> + </detectors> + + <readouts> + <readout name="BubbleDeviceHits"> + <id>system:6,x:12:-6,y:24:-6</id> + </readout> + </readouts> + +</lccdd> diff --git a/examples/OpticalSurfaces/compact/ReadMaterialProperties.xml b/examples/OpticalSurfaces/compact/ReadMaterialProperties.xml new file mode 100644 index 000000000..19d91f024 --- /dev/null +++ b/examples/OpticalSurfaces/compact/ReadMaterialProperties.xml @@ -0,0 +1,233 @@ +<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="TestMaterialProperties" + title="Test reading of TGeo's Material Properties" + author="Markus Frank" + url="None" + status="development" + version="1.0"> + <comment>Test reading of TGeo's Material Properties</comment> + </info> + + <define> + <constant name="world_side" value="1*m"/> + <constant name="world_x" value="world_side/2"/> + <constant name="world_y" value="world_side/2"/> + <constant name="world_z" value="world_side/2"/> + </define> + + <properties> + <matrix name="RINDEX__0x123aff00" coldim="2" values=" + 2.034*eV 1.3435 + 2.068*eV 1.344 + 2.103*eV 1.3445 + 2.139*eV 1.345 + 2.177*eV 1.3455 + 2.216*eV 1.346 + 2.256*eV 1.3465 + 2.298*eV 1.347 + 2.341*eV 1.3475 + 2.386*eV 1.348 + 2.433*eV 1.3485 + 2.481*eV 1.3492 + 2.532*eV 1.35 + 2.585*eV 1.3505 + 2.640*eV 1.351 + 2.697*eV 1.3518 + 2.757*eV 1.3522 + 2.820*eV 1.3530 + 2.885*eV 1.3535 + 2.954*eV 1.354 + 3.026*eV 1.3545 + 3.102*eV 1.355 + 3.181*eV 1.3555 + 3.265*eV 1.356 + 3.353*eV 1.3568 + 3.446*eV 1.3572 + 3.545*eV 1.358 + 3.649*eV 1.3585 + 3.760*eV 1.359 + 3.877*eV 1.3595 + 4.002*eV 1.36 + 4.136*eV 1.3608 + "/> + <matrix name="ABSLENGTH__0x123aff00" coldim="2" values=" + 2.034*eV 3.448*m + 2.068*eV 4.082*m + 2.103*eV 6.329*m + 2.139*eV 9.174*m + 2.177*eV 12.346*m + 2.216*eV 13.889*m + 2.256*eV 15.152*m + 2.298*eV 17.241*m + 2.341*eV 18.868*m + 2.386*eV 20.000*m + 2.433*eV 26.316*m + 2.481*eV 35.714*m + 2.532*eV 45.455*m + 2.585*eV 47.619*m + 2.640*eV 52.632*m + 2.697*eV 52.632*m + 2.757*eV 55.556*m + 2.820*eV 52.632*m + 2.885*eV 52.632*m + 2.954*eV 47.619*m + 3.026*eV 45.455*m + 3.102*eV 41.667*m + 3.181*eV 37.037*m + 3.265*eV 33.333*m + 3.353*eV 30.000*m + 3.446*eV 28.500*m + 3.545*eV 27.000*m + 3.649*eV 24.500*m + 3.760*eV 22.000*m + 3.877*eV 19.500*m + 4.002*eV 17.500*m + 4.136*eV 14.500*m + "/> + <matrix name= "FASTCOMPONENT__0x123aff00" coldim="2" values=" + 2.034*eV 1 + 2.068*eV 1 + 2.103*eV 1 + 2.139*eV 1 + 2.177*eV 1 + 2.216*eV 1 + 2.256*eV 1 + 2.298*eV 1 + 2.341*eV 1 + 2.386*eV 1 + 2.433*eV 1 + 2.481*eV 1 + 2.532*eV 1 + 2.585*eV 1 + 2.640*eV 1 + 2.697*eV 1 + 2.757*eV 1 + 2.820*eV 1 + 2.885*eV 1 + 2.954*eV 1 + 3.026*eV 1 + 3.102*eV 1 + 3.181*eV 1 + 3.265*eV 1 + 3.353*eV 1 + 3.446*eV 1 + 3.545*eV 1 + 3.649*eV 1 + 3.760*eV 1 + 3.877*eV 1 + 4.002*eV 1 + 4.136*eV 1 + "/> + <matrix name= "SLOWCOMPONENT__0x123aff00" coldim="2" values=" + 2.034*eV 0.01 + 2.068*eV 1 + 2.103*eV 2 + 2.139*eV 3 + 2.177*eV 4 + 2.216*eV 5 + 2.256*eV 6 + 2.298*eV 7 + 2.341*eV 8 + 2.386*eV 9 + 2.433*eV 8 + 2.481*eV 7 + 2.532*eV 6 + 2.585*eV 4 + 2.640*eV 3 + 2.697*eV 2 + 2.757*eV 1 + 2.820*eV 0.01 + 2.885*eV 1 + 2.954*eV 2 + 3.026*eV 3 + 3.102*eV 4 + 3.181*eV 5 + 3.265*eV 6 + 3.353*eV 7 + 3.446*eV 8 + 3.545*eV 9 + 3.649*eV 8 + 3.760*eV 7 + 3.877*eV 6 + 4.002*eV 5 + 4.136*eV 4 + "/> + <ignore> + <matrix name= "REFLECTIVITY__0x123aff00" coldim="2" values=" + 2.034*eV + 2.068*eV + 2.103*eV + 2.139*eV + 2.177*eV + 2.216*eV + 2.256*eV + 2.298*eV + 2.341*eV + 2.386*eV + 2.433*eV + 2.481*eV + 2.532*eV + 2.585*eV + 2.640*eV + 2.697*eV + 2.757*eV + 2.820*eV + 2.885*eV + 2.954*eV + 3.026*eV + 3.102*eV + 3.181*eV + 3.265*eV + 3.353*eV + 3.446*eV + 3.545*eV + 3.649*eV + 3.760*eV + 3.877*eV + 4.002*eV + 4.136*eV + "/> + </ignore> + </properties> + + <materials> + <element Z="1" formula="H" name="H" > + <atom type="A" unit="g/mol" value="1.00794" /> + </element> + <element Z="8" formula="O" name="O" > + <atom type="A" unit="g/mol" value="15.9994" /> + </element> + <element Z="7" formula="N" name="N" > + <atom type="A" unit="g/mol" value="14.0068" /> + </element> + <element Z="18" formula="Ar" name="Ar" > + <atom type="A" unit="g/mol" value="39.9477" /> + </element> + + <material name="Air"> + <D type="density" unit="g/cm3" value="0.0012"/> + <fraction n="0.754" ref="N"/> + <fraction n="0.234" ref="O"/> + <fraction n="0.012" ref="Ar"/> + </material> + <!-- We model vakuum just as very thin air --> + <material name="Vacuum"> + <D type="density" unit="g/cm3" value="0.0000000001" /> + <fraction n="0.754" ref="N"/> + <fraction n="0.234" ref="O"/> + <fraction n="0.012" ref="Ar"/> + </material> + <material name="Water"> + <D type="density" value="1.0" unit="g/cm3"/> + <composite n="2" ref="H"/> + <composite n="1" ref="O"/> + <property name="RINDEX" ref="RINDEX__0x123aff00"/> + <property name="ABSLENGTH" ref="ABSLENGTH__0x123aff00"/> + <property name="FASTCOMPONENT" ref="FASTCOMPONENT__0x123aff00"/> + <property name="SLOWCOMPONENT" ref="SLOWCOMPONENT__0x123aff00"/> + </material> + </materials> +</lccdd> diff --git a/examples/OpticalSurfaces/compact/ReadOpticalSurfaces.xml b/examples/OpticalSurfaces/compact/ReadOpticalSurfaces.xml index 843f49710..12b930419 100644 --- a/examples/OpticalSurfaces/compact/ReadOpticalSurfaces.xml +++ b/examples/OpticalSurfaces/compact/ReadOpticalSurfaces.xml @@ -2,13 +2,13 @@ xmlns:xs="http://www.w3.org/2001/XMLSchema" xs:noNamespaceSchemaLocation="http://www.lcsim.org/schemas/compact/1.0/compact.xsd"> - <info name="TestGDMLMatrices" - title="Test reading of TGeo's GDML matrices" + <info name="TestOpticalSurfaces" + title="Test reading of TGeo's Optical Surfaces" author="Markus Frank" url="None" status="development" version="1.0"> - <comment>Test reading of TGeo's GDML matrices</comment> + <comment>Test reading of TGeo's Optical Surfaces</comment> </info> <includes> @@ -23,9 +23,10 @@ <constant name="world_z" value="world_side/2"/> <constant name="PhotMomWaveConv" value="1243.125*eV"/> </define> + <debug> - <type name="matrix" value="1"/> - <type name="surface" value="1"/> + <type name="matrix" value="0"/> + <type name="surface" value="0"/> </debug> <properties> @@ -48,13 +49,15 @@ </properties> <surfaces> - <opticalsurface finish="0" model="0" name="OpticalSurface#1" type="1" value="0"> + <comment> For the values of "finish", model and type, see TGeoOpticalSurface.h ! + </comment> + <opticalsurface finish="Rough_LUT" model="DAVIS" name="OpticalSurface#1" type="dielectric_LUTDAVIS" value="0"> <property name="REFLECTIVITY" ref="REFLECTIVITY0x123aff00"/> - <property name="EFFICIENCY" ref="EFFICIENCY0x8b77240"/> + <property name="EFFICIENCY" ref="EFFICIENCY0x8b77240"/> </opticalsurface> - <opticalsurface finish="0" model="0" name="OpticalSurface#2" type="1" value="1"> + <opticalsurface finish="ground" model="unified" name="OpticalSurface#2" type="dielectric_dielectric" value="1"> <property name="REFLECTIVITY" ref="REFLECTIVITY0x123aff00"/> - <property name="EFFICIENCY" ref="EFFICIENCY0x8b77240"/> + <property name="EFFICIENCY" ref="EFFICIENCY0x8b77240"/> </opticalsurface> </surfaces> diff --git a/examples/OpticalSurfaces/scripts/OpNovice.py b/examples/OpticalSurfaces/scripts/OpNovice.py new file mode 100644 index 000000000..b2f5b483b --- /dev/null +++ b/examples/OpticalSurfaces/scripts/OpNovice.py @@ -0,0 +1,116 @@ +# +# +import os, sys, time, DDG4 +from DDG4 import OutputLevel as Output +from SystemOfUnits import * +# +# +""" + + dd4hep simulation example setup using the python configuration + + @author M.Frank + @version 1.0 + +""" +def run(): + kernel = DDG4.Kernel() + install_dir = os.environ['DD4hepExamplesINSTALL'] + kernel.loadGeometry("file:"+install_dir+"/examples/OpticalSurfaces/compact/OpNovice.xml") + + DDG4.importConstants(kernel.detectorDescription(),debug=False) + geant4 = DDG4.Geant4(kernel,tracker='Geant4TrackerCombineAction') + geant4.printDetectors() + # Configure UI + if len(sys.argv)>1: + geant4.setupCshUI(macro=sys.argv[1]) + else: + geant4.setupCshUI() + + # Configure field + field = geant4.setupTrackingField(prt=True) + # Configure Event actions + prt = DDG4.EventAction(kernel,'Geant4ParticlePrint/ParticlePrint') + prt.OutputLevel = Output.DEBUG + prt.OutputType = 3 # Print both: table and tree + kernel.eventAction().adopt(prt) + + generator_output_level = Output.INFO + + # Configure G4 geometry setup + seq,act = geant4.addDetectorConstruction("Geant4DetectorGeometryConstruction/ConstructGeo") + act.DebugMaterials = True + act.DebugElements = False + act.DebugVolumes = True + act.DebugShapes = True + act.DebugSurfaces = True + + # Configure I/O + evt_root = geant4.setupROOTOutput('RootOutput','OpNovice_'+time.strftime('%Y-%m-%d_%H-%M')) + + # Setup particle gun + gun = geant4.setupGun("Gun",particle='e-',energy=2*GeV,multiplicity=1) + gun.OutputLevel = generator_output_level + + # And handle the simulation particles. + part = DDG4.GeneratorAction(kernel,"Geant4ParticleHandler/ParticleHandler") + kernel.generatorAction().adopt(part) + part.SaveProcesses = ['Decay'] + part.MinimalKineticEnergy = 100*MeV + part.OutputLevel = Output.INFO #generator_output_level + part.enableUI() + user = DDG4.Action(kernel,"Geant4TCUserParticleHandler/UserParticleHandler") + user.TrackingVolume_Zmax = 3.0*m + user.TrackingVolume_Rmax = 3.0*m + user.enableUI() + part.adopt(user) + + geant4.setupTracker('BubbleDevice') + + # Now build the physics list: + phys = geant4.setupPhysics('') + ph = DDG4.PhysicsList(kernel,'Geant4OpticalPhotonPhysics/OpticalGammaPhys') + ph.VerboseLevel = 2 + ph.addParticleGroup('G4BosonConstructor') + ph.addParticleGroup('G4LeptonConstructor') + ph.addParticleGroup('G4MesonConstructor') + ph.addParticleGroup('G4BaryonConstructor') + ph.addParticleGroup('G4IonConstructor') + ph.addParticleConstructor('G4OpticalPhoton') + + ph.addDiscreteParticleProcess('gamma', 'G4GammaConversion') + ph.addDiscreteParticleProcess('gamma', 'G4ComptonScattering') + ph.addDiscreteParticleProcess('gamma', 'G4PhotoElectricEffect') + ph.addParticleProcess('e[+-]', 'G4eMultipleScattering', -1, 1, 1) + ph.addParticleProcess('e[+-]', 'G4eIonisation', -1, 2, 2) + ph.addParticleProcess('e[+-]', 'G4eBremsstrahlung', -1, 3, 3) + ph.addParticleProcess('e+', 'G4eplusAnnihilation', 0,-1, 4) + ph.addParticleProcess('mu[+-]', 'G4MuMultipleScattering',-1, 1, 1) + ph.addParticleProcess('mu[+-]', 'G4MuIonisation', -1, 2, 2) + ph.addParticleProcess('mu[+-]', 'G4MuBremsstrahlung', -1, 3, 3) + ph.addParticleProcess('mu[+-]', 'G4MuPairProduction', -1, 4, 4) + ph.enableUI() + phys.adopt(ph) + + ph = DDG4.PhysicsList(kernel,'Geant4ScintillationPhysics/ScintillatorPhys') + ph.ScintillationYieldFactor = 1.0 + ph.ScintillationExcitationRatio = 1.0 + ph.TrackSecondariesFirst = False + ph.VerboseLevel = 2 + ph.enableUI() + phys.adopt(ph) + + ph = DDG4.PhysicsList(kernel,'Geant4CerenkovPhysics/CerenkovPhys') + ph.MaxNumPhotonsPerStep = 10 + ph.MaxBetaChangePerStep = 10.0 + ph.TrackSecondariesFirst = True + ph.VerboseLevel = 2 + ph.enableUI() + phys.adopt(ph) + + phys.dump() + + geant4.execute() + +if __name__ == "__main__": + run() diff --git a/examples/OpticalSurfaces/src/OpNovice_geo.cpp b/examples/OpticalSurfaces/src/OpNovice_geo.cpp new file mode 100644 index 000000000..6d8921758 --- /dev/null +++ b/examples/OpticalSurfaces/src/OpNovice_geo.cpp @@ -0,0 +1,58 @@ +//========================================================================== +// AIDA Detector description implementation +//-------------------------------------------------------------------------- +// 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 +// +//========================================================================== + +// Include files +#include "DD4hep/DetFactoryHelper.h" +#include "DD4hep/OpticalSurfaces.h" +#include "DD4hep/Printout.h" +#include "DD4hep/Detector.h" + +#include <iostream> +#include <map> + +using namespace std; +using namespace dd4hep; +using namespace dd4hep::detail; + +/// Example of a water box with a box of air inside, the "bubble" +static Ref_t create_detector(Detector &description, xml_h e, SensitiveDetector /* sens */) { + xml_det_t x_det = e; + string det_name = x_det.nameStr(); + DetElement sdet(det_name,x_det.id()); + xml_det_t x_tank = x_det.child(_Unicode(tank)); + xml_det_t x_bubble = x_det.child(_Unicode(bubble)); + + Box tank_box(x_tank.x(), x_tank.y(), x_tank.z()); + Volume tank_vol("Tank",tank_box,description.material(x_tank.attr<string>(_U(material)))); + Box bubble_box(x_bubble.x(), x_bubble.y(), x_bubble.z()); + Volume bubble_vol("Bubble",bubble_box,description.material(x_bubble.attr<string>(_U(material)))); + + tank_vol.setVisAttributes(description, x_tank.attr<string>(_U(vis))); + bubble_vol.setVisAttributes(description, x_bubble.attr<string>(_U(vis))); + + PlacedVolume bubblePlace = tank_vol.placeVolume(bubble_vol); + PlacedVolume tankPlace = description.pickMotherVolume(sdet).placeVolume(tank_vol); + sdet.setPlacement(tankPlace); + + // Now attach the surface + OpticalSurfaceManager surfMgr = description.surfaceManager(); + PlacedVolume hallPlace(description.manager().GetTopNode()); + OpticalSurface waterSurf = surfMgr.opticalSurface("/world/"+det_name+"#WaterSurface"); + OpticalSurface airSurf = surfMgr.opticalSurface("/world/"+det_name+"#AirSurface"); + BorderSurface tankSurf = BorderSurface(description, sdet, "HallTank", waterSurf, tankPlace, hallPlace); + BorderSurface bubbleSurf = BorderSurface(description, sdet, "TankBubble", airSurf, bubblePlace, tankPlace); + bubbleSurf.isValid(); + tankSurf.isValid(); + return sdet; +} +DECLARE_DETELEMENT(DD4hep_OpNovice,create_detector) -- GitLab