Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2026-09-19 08:32:13 

0001 //==========================================================================
0002 //  AIDA Detector description implementation 
0003 //--------------------------------------------------------------------------
0004 // Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
0005 // All rights reserved.
0006 //
0007 // For the licensing terms see $DD4hepINSTALL/LICENSE.
0008 // For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
0009 //
0010 // Author     : M.Frank
0011 //
0012 //==========================================================================
0013 
0014 // Framework include files
0015 #include <DD4hep/Shapes.h>
0016 #include <DD4hep/Volumes.h>
0017 #include <DD4hep/Plugins.h>
0018 #include <DD4hep/Printout.h>
0019 #include <DD4hep/Detector.h>
0020 #include <DD4hep/DD4hepUnits.h>
0021 #include <DD4hep/PropertyTable.h>
0022 #include <DD4hep/DetectorTools.h>
0023 #include <DD4hep/detail/ShapesInterna.h>
0024 #include <DD4hep/detail/ObjectsInterna.h>
0025 #include <DD4hep/detail/DetectorInterna.h>
0026 
0027 #include <DDG4/Geant4Field.h>
0028 #include <DDG4/Geant4Helpers.h>
0029 #include <DDG4/Geant4Converter.h>
0030 #include <DDG4/Geant4UserLimits.h>
0031 #include <DDG4/Geant4AssemblyVolume.h>
0032 #include <DDG4/Geant4PlacementParameterisation.h>
0033 #include "Geant4ShapeConverter.h"
0034 
0035 // ROOT includes
0036 #include <TClass.h>
0037 #include <TTimeStamp.h>
0038 #include <TGeoBoolNode.h>
0039 
0040 // Geant4 include files
0041 #include <G4Version.hh>
0042 #include <G4VisAttributes.hh>
0043 #include <G4PVParameterised.hh>
0044 #include <G4ProductionCuts.hh>
0045 #include <G4VUserRegionInformation.hh>
0046 
0047 #include <G4Box.hh>
0048 #include <G4Tubs.hh>
0049 #include <G4Ellipsoid.hh>
0050 #include <G4UnionSolid.hh>
0051 #include <G4ReflectedSolid.hh>
0052 #include <G4SubtractionSolid.hh>
0053 #include <G4IntersectionSolid.hh>
0054 #include <G4VSensitiveDetector.hh>
0055 
0056 #include <G4Region.hh>
0057 #include <G4Element.hh>
0058 #include <G4Isotope.hh>
0059 #include <G4Material.hh>
0060 #include <G4UserLimits.hh>
0061 #include <G4RegionStore.hh>
0062 #include <G4FieldManager.hh>
0063 #include <G4LogicalVolume.hh>
0064 #include <G4OpticalSurface.hh>
0065 #include <G4ReflectionFactory.hh>
0066 #include <G4LogicalSkinSurface.hh>
0067 #include <G4ElectroMagneticField.hh>
0068 #include <G4LogicalBorderSurface.hh>
0069 #include <G4MaterialPropertiesTable.hh>
0070 #if G4VERSION_NUMBER >= 1040
0071 #include <G4MaterialPropertiesIndex.hh>
0072 #endif
0073 #include <G4ScaledSolid.hh>
0074 #include <CLHEP/Units/SystemOfUnits.h>
0075 
0076 // C/C++ include files
0077 #include <iostream>
0078 #include <iomanip>
0079 #include <sstream>
0080 #include <limits>
0081 
0082 namespace units = dd4hep;
0083 using namespace dd4hep::sim;
0084 using namespace dd4hep;
0085 
0086 namespace {
0087 
0088   static constexpr const double CM_2_MM = (CLHEP::centimeter/dd4hep::centimeter);
0089   static constexpr const char* GEANT4_TAG_CUSTOM = "Geant4-custom";
0090   static constexpr const char* GEANT4_TAG_IGNORE = "Geant4-ignore";
0091   static constexpr const char* GEANT4_TAG_PLUGIN = "Geant4-plugin";
0092   static constexpr const char* GEANT4_TAG_BIRKSCONSTANT    = "BirksConstant";
0093   static constexpr const char* GEANT4_TAG_MEE              = "MeanExcitationEnergy";
0094   static constexpr const char* GEANT4_TAG_ENE_PER_ION_PAIR = "MeanEnergyPerIonPair";
0095 
0096   static std::string indent = "";
0097 
0098   template <typename O, typename C, typename F> void handleRefs(const O* o, const C& c, F pmf) {
0099     for (typename C::const_iterator i = c.begin(); i != c.end(); ++i) {
0100       //(o->*pmf)((*i)->GetName(), *i);
0101       (o->*pmf)("", *i);
0102     }
0103   }
0104 
0105   template <typename O, typename C, typename F> void handle(const O* o, const C& c, F pmf) {
0106     for (typename C::const_iterator i = c.begin(); i != c.end(); ++i) {
0107       (o->*pmf)((*i)->GetName(), *i);
0108     }
0109   }
0110 
0111   template <typename O, typename F> void handleArray(const O* o, const TObjArray* c, F pmf) {
0112     TObjArrayIter arr(c);
0113     for(TObject* i = arr.Next(); i; i=arr.Next())
0114       (o->*pmf)(i);
0115   }
0116 
0117   template <typename O, typename C, typename F> void handleMap(const O* o, const C& c, F pmf) {
0118     for (typename C::const_iterator i = c.begin(); i != c.end(); ++i)
0119       (o->*pmf)((*i).first, (*i).second);
0120   }
0121 
0122   template <typename O, typename C, typename F> void handleRMap(const O* o, const C& c, F pmf) {
0123     for (typename C::const_reverse_iterator i = c.rbegin(); i != c.rend(); ++i)  {
0124       //cout << "Handle RMAP [ " << (*i).first << " ]" << std::endl;
0125       handle(o, i->second, pmf);
0126     }
0127   }
0128   template <typename O, typename C, typename F> void handleRMap_(const O* o, const C& c, F pmf) {
0129     for (typename C::const_iterator i = c.begin(); i != c.end(); ++i)  {
0130       const auto& cc = (*i).second;
0131       for (const auto& j : cc)   {
0132         (o->*pmf)(j);
0133       }
0134     }
0135   }
0136 
0137   std::string make_NCName(const std::string& in)   {
0138     std::string res = detail::str_replace(in, "/", "_");
0139     res = detail::str_replace(res, "#", "_");
0140     return res;
0141   }
0142 
0143   bool is_left_handed(const TGeoMatrix* m)   {
0144     const Double_t* r = m->GetRotationMatrix();
0145     if ( r )    {
0146       Double_t det =
0147         r[0]*r[4]*r[8] + r[3]*r[7]*r[2] + r[6]*r[1]*r[5] -
0148         r[2]*r[4]*r[6] - r[5]*r[7]*r[0] - r[8]*r[1]*r[3];
0149       return det < 0e0;
0150     }
0151     return false;
0152   }
0153 
0154   class G4UserRegionInformation : public G4VUserRegionInformation {
0155   public:
0156     Region region;
0157     double threshold;
0158     bool   storeSecondaries;
0159     G4UserRegionInformation()
0160       : threshold(0.0), storeSecondaries(false) {
0161     }
0162     virtual ~G4UserRegionInformation() {
0163     }
0164     virtual void Print() const {
0165       if (region.isValid())
0166         printout(DEBUG, "Region", "Name:%s", region.name());
0167     }
0168   };
0169 
0170   std::pair<double,double> g4PropertyConversion(int index)   {
0171 #if G4VERSION_NUMBER >= 1040
0172     switch(index)  {
0173     case kRINDEX:                         return std::make_pair(CLHEP::keV/units::keV, 1.0);
0174     case kREFLECTIVITY:                   return std::make_pair(CLHEP::keV/units::keV, 1.0);
0175     case kREALRINDEX:                     return std::make_pair(CLHEP::keV/units::keV, 1.0);
0176     case kIMAGINARYRINDEX:                return std::make_pair(CLHEP::keV/units::keV, 1.0);
0177     case kEFFICIENCY:                     return std::make_pair(CLHEP::keV/units::keV, 1.0);
0178     case kTRANSMITTANCE:                  return std::make_pair(CLHEP::keV/units::keV, 1.0);
0179     case kSPECULARLOBECONSTANT:           return std::make_pair(CLHEP::keV/units::keV, 1.0);
0180     case kSPECULARSPIKECONSTANT:          return std::make_pair(CLHEP::keV/units::keV, 1.0);
0181     case kBACKSCATTERCONSTANT:            return std::make_pair(CLHEP::keV/units::keV, 1.0);
0182     case kGROUPVEL:                       return std::make_pair(CLHEP::keV/units::keV, (CLHEP::m/CLHEP::s)/(units::m/units::s));  // meter/second
0183     case kMIEHG:                          return std::make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
0184     case kRAYLEIGH:                       return std::make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);  // ??? says its a length
0185     case kWLSCOMPONENT:                   return std::make_pair(CLHEP::keV/units::keV, 1.0);
0186     case kWLSABSLENGTH:                   return std::make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
0187     case kABSLENGTH:                      return std::make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
0188 #if G4VERSION_NUMBER >= 1100
0189     case kWLSCOMPONENT2:                  return std::make_pair(CLHEP::keV/units::keV, 1.0);
0190     case kWLSABSLENGTH2:                  return std::make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
0191     case kSCINTILLATIONCOMPONENT1:        return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
0192     case kSCINTILLATIONCOMPONENT2:        return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
0193     case kSCINTILLATIONCOMPONENT3:        return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
0194 #else
0195     case kFASTCOMPONENT:                  return std::make_pair(CLHEP::keV/units::keV, 1.0);
0196     case kSLOWCOMPONENT:                  return std::make_pair(CLHEP::keV/units::keV, 1.0);
0197 #endif
0198     case kPROTONSCINTILLATIONYIELD:       return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV); // Yields: 1/energy
0199     case kDEUTERONSCINTILLATIONYIELD:     return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
0200     case kTRITONSCINTILLATIONYIELD:       return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
0201     case kALPHASCINTILLATIONYIELD:        return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
0202     case kIONSCINTILLATIONYIELD:          return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
0203     case kELECTRONSCINTILLATIONYIELD:     return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
0204     default:
0205       break;
0206     }
0207     printout(FATAL,"Geant4Converter", "+++ Cannot convert material property with index: %d", index);
0208 #else
0209     printout(FATAL,"Geant4Converter", "+++ Cannot convert material property with index: %d [Need Geant4 > 10.03]", index);
0210 #endif
0211     return std::make_pair(0e0,0e0);
0212   }
0213 
0214   double g4ConstPropertyConversion(int index)   {
0215 #if G4VERSION_NUMBER >= 1040
0216     switch(index)   {
0217     case kSURFACEROUGHNESS:            return CLHEP::m/units::m;                             // Length
0218     case kISOTHERMAL_COMPRESSIBILITY:  return (CLHEP::m3/CLHEP::keV)/(units::m3/CLHEP::keV); // Volume/Energy
0219     case kRS_SCALE_FACTOR:             return 1.0;  // ??
0220     case kWLSMEANNUMBERPHOTONS:        return 1.0;  // ??
0221     case kWLSTIMECONSTANT:             return CLHEP::second/units::second;                   // Time
0222     case kMIEHG_FORWARD:               return 1.0;
0223     case kMIEHG_BACKWARD:              return 1.0;
0224     case kMIEHG_FORWARD_RATIO:         return 1.0;
0225     case kSCINTILLATIONYIELD:          return units::keV/CLHEP::keV;                         // Energy
0226     case kRESOLUTIONSCALE:             return 1.0;
0227     case kFERMIPOT:                    return CLHEP::keV/units::keV;                         // Energy
0228     case kDIFFUSION:                   return 1.0;
0229     case kSPINFLIP:                    return 1.0;
0230     case kLOSS:                        return 1.0;  // ??
0231     case kLOSSCS:                      return CLHEP::barn/units::barn;  // ??
0232     case kABSCS:                       return CLHEP::barn/units::barn;  // ??
0233     case kSCATCS:                      return CLHEP::barn/units::barn;  // ??
0234     case kMR_NBTHETA:                  return 1.0;
0235     case kMR_NBE:                      return 1.0;
0236     case kMR_RRMS:                     return 1.0;  // ??
0237     case kMR_CORRLEN:                  return CLHEP::m/units::m;                             // Length
0238     case kMR_THETAMIN:                 return 1.0;
0239     case kMR_THETAMAX:                 return 1.0;
0240     case kMR_EMIN:                     return CLHEP::keV/units::keV;                         // Energy
0241     case kMR_EMAX:                     return CLHEP::keV/units::keV;                         // Energy
0242     case kMR_ANGNOTHETA:               return 1.0;
0243     case kMR_ANGNOPHI:                 return 1.0;
0244     case kMR_ANGCUT:                   return 1.0;
0245 
0246 #if G4VERSION_NUMBER >= 1100
0247     case kSCINTILLATIONTIMECONSTANT1:  return CLHEP::second/units::second;                   // Time
0248     case kSCINTILLATIONTIMECONSTANT2:  return CLHEP::second/units::second;                   // Time
0249     case kSCINTILLATIONTIMECONSTANT3:  return CLHEP::second/units::second;                   // Time
0250     case kSCINTILLATIONRISETIME1:      return CLHEP::second/units::second;                   // Time
0251     case kSCINTILLATIONRISETIME2:      return CLHEP::second/units::second;                   // Time
0252     case kSCINTILLATIONRISETIME3:      return CLHEP::second/units::second;                   // Time
0253     case kSCINTILLATIONYIELD1:         return 1.0;
0254     case kSCINTILLATIONYIELD2:         return 1.0;
0255     case kSCINTILLATIONYIELD3:         return 1.0;
0256     case kPROTONSCINTILLATIONYIELD1:   return 1.0;
0257     case kPROTONSCINTILLATIONYIELD2:   return 1.0;
0258     case kPROTONSCINTILLATIONYIELD3:   return 1.0;
0259     case kDEUTERONSCINTILLATIONYIELD1: return 1.0;
0260     case kDEUTERONSCINTILLATIONYIELD2: return 1.0;
0261     case kDEUTERONSCINTILLATIONYIELD3: return 1.0;
0262     case kALPHASCINTILLATIONYIELD1:    return 1.0;
0263     case kALPHASCINTILLATIONYIELD2:    return 1.0;
0264     case kALPHASCINTILLATIONYIELD3:    return 1.0;
0265     case kIONSCINTILLATIONYIELD1:      return 1.0;
0266     case kIONSCINTILLATIONYIELD2:      return 1.0;
0267     case kIONSCINTILLATIONYIELD3:      return 1.0;
0268     case kELECTRONSCINTILLATIONYIELD1: return 1.0;
0269     case kELECTRONSCINTILLATIONYIELD2: return 1.0;
0270     case kELECTRONSCINTILLATIONYIELD3: return 1.0;
0271 #else
0272     case kFASTTIMECONSTANT:            return CLHEP::second/units::second;                   // Time
0273     case kFASTSCINTILLATIONRISETIME:   return CLHEP::second/units::second;                   // Time
0274     case kSLOWTIMECONSTANT:            return CLHEP::second/units::second;                   // Time
0275     case kSLOWSCINTILLATIONRISETIME:   return CLHEP::second/units::second;                   // Time
0276     case kYIELDRATIO:                  return 1.0;
0277 #endif
0278     default:
0279       break;
0280     }
0281     printout(FATAL,"Geant4Converter", "+++ Cannot convert CONST material property with index: %d", index);
0282 #else
0283     printout(FATAL,"Geant4Converter", "+++ Cannot convert material property with index: %d [Need Geant4 > 10.03]", index);
0284 #endif
0285     return 0.0;
0286   }
0287 }
0288 
0289 /// Initializing Constructor
0290 Geant4Converter::Geant4Converter(const Detector& description_ref)
0291   : Geant4Mapping(description_ref), checkOverlaps(true) {
0292   this->Geant4Mapping::init();
0293   m_propagateRegions = true;
0294   outputLevel = PrintLevel(printLevel() - 1);
0295 }
0296 
0297 /// Initializing Constructor
0298 Geant4Converter::Geant4Converter(const Detector& description_ref, PrintLevel level)
0299   : Geant4Mapping(description_ref), outputLevel(level)  {
0300   this->Geant4Mapping::init();
0301   m_propagateRegions = true;
0302 }
0303 
0304 /// Standard destructor
0305 Geant4Converter::~Geant4Converter() {
0306 }
0307 
0308 /// Handle the conversion of isotopes
0309 void* Geant4Converter::handleIsotope(const std::string& /* name */, const TGeoIsotope* iso) const {
0310   G4Isotope* g4i = data().g4Isotopes[iso];
0311   if ( !g4i )  {
0312     double a_conv = (CLHEP::g / CLHEP::mole);
0313     g4i = new G4Isotope(iso->GetName(), iso->GetZ(), iso->GetN(), iso->GetA()*a_conv);
0314     printout(debugElements ? ALWAYS : outputLevel,
0315              "Geant4Converter", "++ Created G4 Isotope %s from data: Z=%d N=%d A=%.3f [g/mole]",
0316              iso->GetName(), iso->GetZ(), iso->GetN(), iso->GetA());
0317     data().g4Isotopes[iso] = g4i;
0318   }
0319   return g4i;
0320 }
0321 
0322 /// Handle the conversion of elements
0323 void* Geant4Converter::handleElement(const std::string& name, const Atom element) const {
0324   G4Element* g4e = data().g4Elements[element];
0325   if ( !g4e ) {
0326     PrintLevel lvl = debugElements ? ALWAYS : outputLevel;
0327     if (element->GetNisotopes() > 0) {
0328       g4e = new G4Element(name, element->GetTitle(), element->GetNisotopes());
0329       for (int i = 0, n = element->GetNisotopes(); i < n; ++i) {
0330         TGeoIsotope* iso = element->GetIsotope(i);
0331         G4Isotope* g4iso = (G4Isotope*)handleIsotope(iso->GetName(), iso);
0332         g4e->AddIsotope(g4iso, element->GetRelativeAbundance(i));
0333       }
0334     }
0335     else {
0336       // This adds in Geant4 the natural isotopes, which we normally do not want. We want to steer it outselves.
0337       double a_conv = (CLHEP::g / CLHEP::mole);
0338       g4e = new G4Element(element->GetTitle(), name, element->Z(), element->A()*a_conv);
0339       printout(lvl, "Geant4Converter", "++ Created G4 Isotope %s from data: Z=%d N=%d A=%.3f [g/mole]",
0340                element->GetName(), element->Z(), element->N(), element->A());
0341     }
0342     std::stringstream str;
0343     str << (*g4e) << std::endl;
0344     printout(lvl, "Geant4Converter", "++ Created G4 element %s", str.str().c_str());
0345     data().g4Elements[element] = g4e;
0346   }
0347   return g4e;
0348 }
0349 
0350 /// Dump material in GDML format to output stream
0351 void* Geant4Converter::handleMaterial(const std::string& name, Material medium) const {
0352   Geant4GeometryInfo& info = data();
0353   G4Material*         mat  = info.g4Materials[medium];
0354   if ( !mat )  {
0355     PrintLevel    lvl      = debugMaterials ? ALWAYS : outputLevel;
0356     TGeoMaterial* material = medium->GetMaterial();
0357     G4State       state    = kStateUndefined;
0358     double        density  = material->GetDensity() * (CLHEP::gram / CLHEP::cm3);
0359     if ( density < 1e-25 )
0360       density = 1e-25;
0361     switch ( material->GetState() ) {
0362     case TGeoMaterial::kMatStateSolid:
0363       state = kStateSolid;
0364       break;
0365     case TGeoMaterial::kMatStateLiquid:
0366       state = kStateLiquid;
0367       break;
0368     case TGeoMaterial::kMatStateGas:
0369       state = kStateGas;
0370       break;
0371     default:
0372     case TGeoMaterial::kMatStateUndefined:
0373       state = kStateUndefined;
0374       break;
0375     }
0376     printout(lvl,"Geant4Material","+++ Setting up material %s", name.c_str());
0377     if ( material->IsMixture() )  {
0378       double A_total = 0.0;
0379       double W_total = 0.0;
0380       TGeoMixture* mix = (TGeoMixture*) material;
0381       int    nElements = mix->GetNelements();
0382       mat = new G4Material(name, density, nElements, state, 
0383                            material->GetTemperature(), material->GetPressure());
0384       for (int i = 0; i < nElements; ++i)  {
0385         A_total += (mix->GetAmixt())[i];
0386         W_total += (mix->GetWmixt())[i];
0387       }
0388       for (int i = 0; i < nElements; ++i) {
0389         TGeoElement* e = mix->GetElement(i);
0390         G4Element* g4e = (G4Element*) handleElement(e->GetName(), Atom(e));
0391         if (!g4e) {
0392           printout(ERROR, name, 
0393                    "Missing element component %s for material %s. A=%f W=%f", 
0394                    e->GetName(), mix->GetName(), A_total, W_total);
0395         }
0396         //mat->AddElement(g4e, (mix->GetAmixt())[i] / A_total);
0397         mat->AddElement(g4e, (mix->GetWmixt())[i] / W_total);
0398       }
0399     }
0400     else {
0401       double z = material->GetZ(), a = material->GetA();
0402       if ( z < 1.0000001 ) z = 1.0;
0403       if ( a < 0.5000001 ) a = 1.0;
0404       mat = new G4Material(name, z, a, density, state, 
0405                            material->GetTemperature(), material->GetPressure());
0406     }
0407 
0408     std::string plugin_name { };
0409     double value = 0e0;
0410     double ionisation_mee = -2e100;
0411     double ionisation_birks_constant = -2e100;
0412     double ionisation_ene_per_ion_pair = -2e100;
0413 
0414     /// Attach the material properties if any
0415     G4MaterialPropertiesTable* tab = 0;
0416     TListIter propIt(&material->GetProperties());
0417     for(TObject* obj=propIt.Next(); obj; obj = propIt.Next())  {
0418       std::string  exc_str;
0419       bool         custom_property = false;
0420       TNamed*      named  = (TNamed*)obj;
0421       TGDMLMatrix* matrix = info.manager->GetGDMLMatrix(named->GetTitle());
0422       const char*  cptr   = ::strstr(matrix->GetName(), GEANT4_TAG_IGNORE);
0423       if( nullptr != cptr )  {
0424         printout(INFO,name,"++ Ignore property %s [%s]. Not Suitable for Geant4.",
0425                  matrix->GetName(), matrix->GetTitle());
0426         continue;
0427       }
0428       cptr = ::strstr(matrix->GetTitle(), GEANT4_TAG_IGNORE);
0429       if( nullptr != cptr )  {
0430         printout(INFO,name,"++ Ignore property %s [%s]. Not Suitable for Geant4.",
0431                  matrix->GetName(), matrix->GetTitle());
0432         continue;
0433       }
0434       cptr   = ::strstr(matrix->GetName(), GEANT4_TAG_CUSTOM);
0435       if( nullptr != cptr )  {
0436         custom_property = true;
0437       }
0438       cptr = ::strstr(matrix->GetTitle(), GEANT4_TAG_CUSTOM);
0439       if( nullptr != cptr )  {
0440         custom_property = true;
0441       }
0442 
0443       Geant4GeometryInfo::PropertyVector* v =
0444         (Geant4GeometryInfo::PropertyVector*)handleMaterialProperties(matrix);
0445       if( nullptr == v )  {
0446         except("Geant4Converter", "++ FAILED to create G4 material %s [Cannot convert property:%s]",
0447                material->GetName(), named->GetName());
0448       }
0449       if( nullptr == tab )  {
0450         tab = new G4MaterialPropertiesTable();
0451         mat->SetMaterialPropertiesTable(tab);
0452       }
0453       int idx = -1;
0454       try  {
0455         if( !custom_property )  {
0456           const auto& pn = tab->GetMaterialPropertyNames();
0457           if( std::find(std::begin(pn), std::end(pn), named->GetName()) != pn.end() )  {
0458             idx = tab->GetPropertyIndex(named->GetName());
0459           }
0460         }
0461       }
0462       catch(const std::exception& e)   {
0463         exc_str = e.what();
0464         idx = -1;
0465       }
0466       catch(...)   {
0467         idx = -1;
0468       }
0469       if ( idx < 0 && !custom_property )  {
0470         printout(ERROR, "Geant4Converter",
0471                  "++ UNKNOWN Geant4 Property: %-20s %s [IGNORED]",
0472                  exc_str.c_str(), named->GetName());
0473         continue;
0474       }
0475       // We need to convert the property from TGeo units to Geant4 units
0476       std::vector<double> bins(v->bins), vals(v->values);
0477       std::pair<double, double> conv = { 1e0, 1e0 };
0478       if( !custom_property )  {
0479         conv = g4PropertyConversion(idx);
0480         for(std::size_t i=0, count=bins.size(); i<count; ++i)
0481           bins[i] *= conv.first, vals[i] *= conv.second;
0482       }
0483       G4MaterialPropertyVector* vec =
0484         new G4MaterialPropertyVector(&bins[0], &vals[0], bins.size());
0485 #if G4VERSION_NUMBER >= 1100
0486       tab->AddProperty(named->GetName(), vec, custom_property);
0487 #else
0488       tab->AddProperty(named->GetName(), vec);
0489 #endif
0490       printout(lvl, name, "++      %sProperty: %-20s [%ld x %ld] -> %s ",
0491                custom_property ? "CUSTOM " : "", named->GetName(),
0492                matrix->GetRows(), matrix->GetCols(), named->GetTitle());
0493       for(std::size_t i=0, count=v->bins.size(); i<count; ++i)
0494         printout(lvl, name, "  Geant4: %s %8.3g [MeV]  TGeo: %8.3g [GeV] Conversion: %8.3g",
0495                  named->GetName(), bins[i], v->bins[i], conv.first);
0496     }
0497 
0498     /// Attach the material properties if any
0499     TListIter cpropIt(&material->GetConstProperties());
0500     for(TObject* obj=cpropIt.Next(); obj; obj = cpropIt.Next())  {
0501       std::string  exc_str;
0502       Bool_t  err = kFALSE;
0503       TNamed* named = (TNamed*)obj;
0504       bool    custom_property = false;
0505 
0506       const char*  cptr = ::strstr(named->GetName(), GEANT4_TAG_IGNORE);
0507       if( nullptr != cptr )   {
0508         printout(INFO, name, "++ Ignore CONST property %s [%s].",
0509                  named->GetName(), named->GetTitle());
0510         continue;
0511       }
0512       cptr = ::strstr(named->GetTitle(), GEANT4_TAG_IGNORE);
0513       if( nullptr != cptr )  {
0514         printout(INFO, name,"++ Ignore CONST property %s [%s].",
0515                  named->GetName(), named->GetTitle());
0516         continue;
0517       }
0518       cptr = ::strstr(named->GetName(), GEANT4_TAG_PLUGIN);
0519       if( nullptr != cptr )  {
0520         printout(INFO, name, "++ Ignore CONST property %s [%s]  --> Plugin.",
0521                  named->GetName(), named->GetTitle());
0522         plugin_name = named->GetTitle();
0523         continue;
0524       }
0525       cptr = ::strstr(named->GetName(), GEANT4_TAG_BIRKSCONSTANT);
0526       if( nullptr != cptr )  {
0527         err = kFALSE;
0528         value = material->GetConstProperty(GEANT4_TAG_BIRKSCONSTANT,&err);
0529         if ( err == kFALSE ) ionisation_birks_constant = value * (CLHEP::mm/CLHEP::MeV)/(units::mm/units::MeV);
0530         continue;
0531       }
0532       cptr = ::strstr(named->GetName(), GEANT4_TAG_MEE);
0533       if( nullptr != cptr )  {
0534         err = kFALSE;
0535         value = material->GetConstProperty(GEANT4_TAG_MEE, &err);
0536         if ( err == kFALSE ) ionisation_mee = value * (CLHEP::MeV/units::MeV);
0537         continue;
0538       }
0539       cptr = ::strstr(named->GetName(), GEANT4_TAG_ENE_PER_ION_PAIR);
0540       if( nullptr != cptr )  {
0541         err = kFALSE;
0542         value = material->GetConstProperty(GEANT4_TAG_ENE_PER_ION_PAIR,&err);
0543         if ( err == kFALSE ) ionisation_ene_per_ion_pair = value * (CLHEP::MeV/units::MeV);
0544         continue;
0545       }
0546       cptr = ::strstr(named->GetName(), GEANT4_TAG_CUSTOM);
0547       if ( nullptr != cptr )  {
0548         custom_property = true;
0549       }
0550       cptr = ::strstr(named->GetTitle(), GEANT4_TAG_CUSTOM);
0551       if ( nullptr != cptr )  {
0552         custom_property = true;
0553       }
0554 
0555       err = kFALSE;
0556       value = info.manager->GetProperty(named->GetTitle(), &err);
0557       if ( err != kFALSE )  {
0558         except(name,
0559                "++ FAILED to create G4 material %s [Cannot convert const property: %s]",
0560                material->GetName(), named->GetName());
0561       }
0562       if ( nullptr == tab )  {
0563         tab = new G4MaterialPropertiesTable();
0564         mat->SetMaterialPropertiesTable(tab);
0565       }
0566       int idx = -1;
0567       try   {
0568         if( !custom_property )  {
0569           const auto& pn = tab->GetMaterialConstPropertyNames();
0570           if( std::find(std::begin(pn), std::end(pn), named->GetName()) != pn.end() )  {
0571             idx = tab->GetConstPropertyIndex(named->GetName());
0572           }
0573         }
0574       }
0575       catch(const std::exception& e)   {
0576         exc_str = e.what();
0577         idx = -1;
0578       }
0579       catch(...)   {
0580         idx = -1;
0581       }
0582       if ( idx < 0 && !custom_property )  {
0583         printout(ERROR, name,
0584                  "++ UNKNOWN Geant4 CONST Property: %-20s %s [IGNORED]",
0585                  exc_str.c_str(), named->GetName());
0586         continue;
0587       }
0588       // We need to convert the property from TGeo units to Geant4 units
0589       if ( !custom_property )  {
0590         double conv = g4ConstPropertyConversion(idx);
0591         value = value * conv;
0592       }
0593       printout(lvl, name, "++      %sCONST Property: %-20s %g ",
0594                custom_property ? "CUSTOM " : "", named->GetName(), value);
0595       tab->AddConstProperty(named->GetName(), value);
0596     }
0597     //
0598     // Set Birk's constant if it was supplied in the material table of the TGeoMaterial
0599     auto* ionisation = mat->GetIonisation();
0600     std::stringstream str;
0601     str << (*mat);
0602     if ( ionisation )   {
0603       if ( ionisation_birks_constant > 0e0 )   {
0604         ionisation->SetBirksConstant(ionisation_birks_constant);
0605       }
0606       if ( ionisation_mee > -1e100 )   {
0607         ionisation->SetMeanExcitationEnergy(ionisation_mee);
0608       }
0609       if ( ionisation_ene_per_ion_pair > 0e0 )   {
0610         ionisation->SetMeanEnergyPerIonPair(ionisation_ene_per_ion_pair);
0611       }
0612       str << "          log(MEE): " << std::setprecision(4) << ionisation->GetLogMeanExcEnergy();
0613       if ( ionisation_birks_constant > 0e0 )
0614         str << "  Birk's constant: " << std::setprecision(4) << ionisation->GetBirksConstant() << " [mm/MeV]";
0615       if ( ionisation_ene_per_ion_pair > 0e0 )
0616         str << "  Mean Energy Per Ion Pair: " << std::setprecision(4) << ionisation->GetMeanEnergyPerIonPair()/CLHEP::eV << " [eV]";
0617     }
0618     else  {
0619       str << "          No ionisation parameters available.";
0620     }
0621     printout(lvl, name, "++ Created G4 material %s", str.str().c_str());
0622 
0623     if ( !plugin_name.empty() )    {
0624       // Call plugin to create extended material if requested
0625       Detector* det = const_cast<Detector*>(&m_detDesc);
0626       G4Material* extended_mat = PluginService::Create<G4Material*>(plugin_name, det, medium, mat);
0627       if ( !extended_mat )   {
0628         except("G4Cnv::material["+name+"]","++ FATAL Failed to call plugin to create material.");
0629       }
0630       mat = extended_mat;
0631     }
0632     info.g4Materials[medium] = mat;
0633   }
0634   return mat;
0635 }
0636 
0637 /// Dump solid in GDML format to output stream
0638 void* Geant4Converter::handleSolid(const std::string& name, const TGeoShape* shape) const {
0639   G4VSolid* solid = nullptr;
0640   if ( shape ) {
0641     if ( nullptr != (solid = data().g4Solids[shape]) )   {
0642       return solid;
0643     }
0644     TClass*    isa = shape->IsA();
0645     PrintLevel lvl = debugShapes ? ALWAYS : outputLevel;
0646     if (isa == TGeoShapeAssembly::Class()) {
0647       // Assemblies have no corresponding 'shape' in Geant4. Ignore the shape translation.
0648       // It does not harm, since this 'shape' is never accessed afterwards.
0649       data().g4Solids[shape] = solid = convertShape<TGeoShapeAssembly>(shape);
0650       return solid;
0651     }
0652     else if (isa == TGeoBBox::Class())
0653       solid = convertShape<TGeoBBox>(shape);
0654     else if (isa == TGeoTube::Class())
0655       solid = convertShape<TGeoTube>(shape);
0656     else if (isa == TGeoTubeSeg::Class())
0657       solid = convertShape<TGeoTubeSeg>(shape);
0658     else if (isa == TGeoCtub::Class())
0659       solid = convertShape<TGeoCtub>(shape);
0660     else if (isa == TGeoEltu::Class())
0661       solid = convertShape<TGeoEltu>(shape);
0662     else if (isa == TwistedTubeObject::Class())
0663       solid = convertShape<TwistedTubeObject>(shape);
0664     else if (isa == TGeoTrd1::Class())
0665       solid = convertShape<TGeoTrd1>(shape);
0666     else if (isa == TGeoTrd2::Class())
0667       solid = convertShape<TGeoTrd2>(shape);
0668     else if (isa == TGeoHype::Class())
0669       solid = convertShape<TGeoHype>(shape);
0670     else if (isa == TGeoXtru::Class())
0671       solid = convertShape<TGeoXtru>(shape);
0672     else if (isa == TGeoPgon::Class())
0673       solid = convertShape<TGeoPgon>(shape);
0674     else if (isa == TGeoPcon::Class())
0675       solid = convertShape<TGeoPcon>(shape);
0676     else if (isa == TGeoCone::Class())
0677       solid = convertShape<TGeoCone>(shape);
0678     else if (isa == TGeoConeSeg::Class())
0679       solid = convertShape<TGeoConeSeg>(shape);
0680     else if (isa == TGeoParaboloid::Class())
0681       solid = convertShape<TGeoParaboloid>(shape);
0682     else if (isa == TGeoSphere::Class())
0683       solid = convertShape<TGeoSphere>(shape);
0684     else if (isa == TGeoTorus::Class())
0685       solid = convertShape<TGeoTorus>(shape);
0686     else if (isa == TGeoTrap::Class())
0687       solid = convertShape<TGeoTrap>(shape);
0688     else if (isa == TGeoArb8::Class()) 
0689       solid = convertShape<TGeoArb8>(shape);
0690     else if (isa == TGeoPara::Class())
0691       solid = convertShape<TGeoPara>(shape);
0692     else if (isa == TGeoTessellated::Class()) 
0693       solid = convertShape<TGeoTessellated>(shape);
0694     else if (isa == TGeoScaledShape::Class())  {
0695       TGeoScaledShape* sh   = (TGeoScaledShape*) shape;
0696       TGeoShape*       sol  = sh->GetShape();
0697       if ( sol->IsA() == TGeoShapeAssembly::Class() )  {
0698         return solid;
0699       }
0700       const double*    vals = sh->GetScale()->GetScale();
0701       G4Scale3D        scal(vals[0], vals[1], vals[2]);
0702       G4VSolid* g4solid = (G4VSolid*)handleSolid(sol->GetName(), sol);
0703       if ( scal.xx()>0e0 && scal.yy()>0e0 && scal.zz()>0e0 )
0704         solid = new G4ScaledSolid(sh->GetName(), g4solid, scal);
0705       else
0706         solid = new G4ReflectedSolid(g4solid->GetName()+"_refl", g4solid, scal);
0707     }
0708     else if ( isa == TGeoCompositeShape::Class() )   {
0709       const TGeoCompositeShape* sh = (const TGeoCompositeShape*) shape;
0710       const TGeoBoolNode* boolean = sh->GetBoolNode();
0711       TGeoBoolNode::EGeoBoolType oper = boolean->GetBooleanOperator();
0712       TGeoMatrix* matrix = boolean->GetRightMatrix();
0713       G4VSolid* left  = (G4VSolid*) handleSolid(name + "_left", boolean->GetLeftShape());
0714       G4VSolid* right = (G4VSolid*) handleSolid(name + "_right", boolean->GetRightShape());
0715       
0716       if (!left) {
0717         except("Geant4Converter","++ No left Geant4 Solid present for composite shape: %s",name.c_str());
0718       }
0719       if (!right) {
0720         except("Geant4Converter","++ No right Geant4 Solid present for composite shape: %s",name.c_str());
0721       }
0722 
0723       TGeoShape* ls = boolean->GetLeftShape();
0724       TGeoShape* rs = boolean->GetRightShape();
0725       if (strcmp(ls->ClassName(), "TGeoScaledShape") == 0 &&
0726           strcmp(rs->ClassName(), "TGeoBBox") == 0) {
0727         if (strcmp(((TGeoScaledShape *)ls)->GetShape()->ClassName(), "TGeoSphere") == 0) {
0728           if (oper == TGeoBoolNode::kGeoIntersection) {
0729             TGeoScaledShape* lls = (TGeoScaledShape *)ls;
0730             TGeoBBox* rrs = (TGeoBBox*)rs;
0731             double sx     = lls->GetScale()->GetScale()[0];
0732             double sy     = lls->GetScale()->GetScale()[1];
0733             double radius = ((TGeoSphere *)lls->GetShape())->GetRmax();
0734             double dz     = rrs->GetDZ();
0735             double zorig  = rrs->GetOrigin()[2];
0736             double zcut2  = dz + zorig;
0737             double zcut1  = 2 * zorig - zcut2;
0738             solid = new G4Ellipsoid(name,
0739                                     sx * radius * CM_2_MM,
0740                                     sy * radius * CM_2_MM,
0741                                     radius * CM_2_MM,
0742                                     zcut1 * CM_2_MM,
0743                                     zcut2 * CM_2_MM);
0744             data().g4Solids[shape] = solid;
0745             return solid;
0746           }
0747         }
0748       }
0749 
0750       if ( matrix->IsRotation() ) {
0751         G4Transform3D transform;
0752         g4Transform(matrix, transform);
0753         if (oper == TGeoBoolNode::kGeoSubtraction)
0754           solid = new G4SubtractionSolid(name, left, right, transform);
0755         else if (oper == TGeoBoolNode::kGeoUnion)
0756           solid = new G4UnionSolid(name, left, right, transform);
0757         else if (oper == TGeoBoolNode::kGeoIntersection)
0758           solid = new G4IntersectionSolid(name, left, right, transform);
0759       }
0760       else {
0761         const Double_t *t = matrix->GetTranslation();
0762         G4ThreeVector transform(t[0] * CM_2_MM, t[1] * CM_2_MM, t[2] * CM_2_MM);
0763         if (oper == TGeoBoolNode::kGeoSubtraction)
0764           solid = new G4SubtractionSolid(name, left, right, 0, transform);
0765         else if (oper == TGeoBoolNode::kGeoUnion)
0766           solid = new G4UnionSolid(name, left, right, 0, transform);
0767         else if (oper == TGeoBoolNode::kGeoIntersection)
0768           solid = new G4IntersectionSolid(name, left, right, 0, transform);
0769       }
0770     }
0771 
0772     if ( !solid )
0773       except("Geant4Converter","++ Failed to handle unknown solid shape: %s of type %s",
0774              name.c_str(), isa->GetName());
0775     printout(lvl,"Geant4Converter","++ Successessfully converted shape [%p] of type:%s to %s.",
0776              solid,isa->GetName(),typeName(typeid(*solid)).c_str());
0777     data().g4Solids[shape] = solid;
0778   }
0779   return solid;
0780 }
0781 
0782 /// Dump logical volume in GDML format to output stream
0783 void* Geant4Converter::handleVolume(const std::string& name, const TGeoVolume* volume) const {
0784   Volume _v(volume);
0785   Geant4GeometryInfo& info = data();
0786   PrintLevel lvl = debugVolumes ? ALWAYS : outputLevel;
0787   Geant4GeometryMaps::VolumeMap::const_iterator volIt = info.g4Volumes.find(volume);
0788   if ( _v.testFlagBit(Volume::VETO_SIMU) )  {
0789     printout(lvl, "Geant4Converter",
0790              "++ Volume %s not converted [Veto'ed for simulation]",
0791              volume->GetName());
0792     return nullptr;
0793   }
0794   else if (volIt == info.g4Volumes.end() ) {
0795     const char*  vnam = volume->GetName();
0796     TGeoMedium*  med  = volume->GetMedium();
0797     Solid        sh   = volume->GetShape();
0798     bool         is_assembly = sh->IsA() == TGeoShapeAssembly::Class() || volume->IsAssembly();
0799 
0800     printout(lvl, "Geant4Converter", "++ Convert Volume %-32s: %p %s/%s assembly:%s",
0801              vnam, volume, sh.type(), _v.type(), yes_no(is_assembly));
0802     if ( is_assembly ) {
0803       return nullptr;
0804     }
0805     Region        reg      = _v.region();
0806     LimitSet      lim      = _v.limitSet();
0807     VisAttr       vis      = _v.visAttributes();
0808     G4Region*     g4region = reg.isValid() ? info.g4Regions[reg] : nullptr;
0809     G4UserLimits* g4limits = lim.isValid() ? info.g4Limits[lim]  : nullptr;
0810     G4VSolid*     g4solid  = (G4VSolid*)   handleSolid(sh->GetName(), sh);
0811     G4Material*   g4medium = (G4Material*) handleMaterial(med->GetName(), Material(med));
0812     /// Check all pre-conditions
0813     if ( !g4solid )   {
0814       except("G4Converter","++ No Geant4 Solid present for volume: %s", vnam);
0815     }
0816     else if ( !g4medium )   {
0817       except("G4Converter","++ No Geant4 material present for volume: %s", vnam);
0818     }
0819     else if ( reg.isValid() && !g4region )  {
0820       except("G4Cnv::volume["+name+"]"," ++ Failed to access Geant4 region %s.", reg.name());
0821     }
0822     else if ( lim.isValid() && !g4limits )  {
0823       except("G4Cnv::volume["+name+"]","++ FATAL Failed to access Geant4 user limits %s.", lim.name());
0824     }
0825     else if ( g4limits )   {
0826       printout(lvl, "Geant4Converter", "++ Volume     + Apply LIMITS settings: %-24s to volume %s.",
0827                lim.name(), vnam);
0828     }
0829 
0830     G4LogicalVolume* g4vol = nullptr;
0831     if( _v.hasProperties() && !_v.getProperty(GEANT4_TAG_PLUGIN,"").empty() )   {
0832       Detector*   det = const_cast<Detector*>(&m_detDesc); 
0833       std::string plugin = _v.getProperty(GEANT4_TAG_PLUGIN,"");
0834       g4vol = PluginService::Create<G4LogicalVolume*>(plugin, det, _v, g4solid, g4medium);
0835       if ( !g4vol )    {
0836         except("G4Cnv::volume["+name+"]","++ FATAL Failed to call plugin to create logical volume.");
0837       }
0838     }
0839     else  {
0840       g4vol = new G4LogicalVolume(g4solid, g4medium, vnam, nullptr, nullptr, nullptr);
0841     }
0842     PrintLevel plevel = (debugVolumes||debugRegions||debugLimits) ? ALWAYS : outputLevel;
0843     /// Set smartless optimization
0844     unsigned char smart_less_value = _v.smartlessValue();
0845     if( smart_less_value != Volume::NO_SMARTLESS_OPTIMIZATION )  {
0846       printout(ALWAYS, "Geant4Converter",
0847                "++ Volume %s Set Smartless value to %d",
0848                vnam, int(smart_less_value));
0849       g4vol->SetSmartless( smart_less_value );
0850     }
0851     /// Assign limits if necessary
0852     if( g4limits )   {
0853       g4vol->SetUserLimits(g4limits);
0854     }
0855     if( g4region )   {
0856       printout(plevel, "Geant4Converter",
0857                "++ Volume     + Apply REGION settings: %-24s to volume %s.",
0858                reg.name(), vnam);
0859       // Handle the region settings for the world volume seperately.
0860       // Geant4 does NOT WANT any regions assigned to the workd volume.
0861       // The world's region is created in the G4RunManagerKernel!
0862       if ( _v == m_detDesc.worldVolume() )   {
0863         const char* wrd_nam = "DefaultRegionForTheWorld";
0864         const char* src_nam = g4region->GetName().c_str();
0865         auto* world_region  = G4RegionStore::GetInstance()->GetRegion(wrd_nam, false);
0866         if ( auto* cuts = g4region->GetProductionCuts() )   {
0867           world_region->SetProductionCuts(cuts);
0868           printout(plevel, "Geant4Converter",
0869                    "++ Volume %s Region: %s. Apply production cuts from %s", 
0870                    vnam, wrd_nam, src_nam);
0871         }
0872         if ( auto* lims = g4region->GetUserLimits() )   {
0873           world_region->SetUserLimits(lims);
0874           printout(plevel, "Geant4Converter",
0875                    "++ Volume %s Region: %s. Apply user limits from %s", 
0876                    vnam, wrd_nam, src_nam);
0877         }
0878       }
0879       else   {
0880         g4vol->SetRegion(g4region);
0881         g4region->AddRootLogicalVolume(g4vol);
0882       }
0883     }
0884     G4VisAttributes* g4vattr = vis.isValid()
0885       ? (G4VisAttributes*)handleVis(vis.name(), vis) : nullptr;
0886     if ( g4vattr )   {
0887       g4vol->SetVisAttributes(g4vattr);
0888     }
0889     info.g4Volumes[volume] = g4vol;
0890     printout(lvl, "Geant4Converter",
0891              "++ Volume     + %s converted: %p ---> G4: %p", vnam, volume, g4vol);
0892   }
0893   return nullptr;
0894 }
0895 
0896 /// Dump logical volume in GDML format to output stream
0897 void* Geant4Converter::collectVolume(const std::string& /* name */, const TGeoVolume* volume) const {
0898   Geant4GeometryInfo& info = data();
0899   Volume              _v(volume);
0900   Region              reg = _v.region();
0901   LimitSet            lim = _v.limitSet();
0902   SensitiveDetector   det = _v.sensitiveDetector();
0903   bool              world = (volume == m_detDesc.worldVolume().ptr());
0904 
0905   if ( !world )   {
0906     if ( lim.isValid() )
0907       info.limits[lim].insert(volume);
0908     if ( reg.isValid() )
0909       info.regions[reg].insert(volume);
0910     if ( det.isValid() )
0911       info.sensitives[det].insert(volume);
0912   }
0913   return (void*)volume;
0914 }
0915 
0916 /// Dump volume placement in GDML format to output stream
0917 void* Geant4Converter::handleAssembly(const std::string& name, const TGeoNode* node) const {
0918   TGeoVolume* mot_vol = node->GetVolume();
0919   PrintLevel lvl = debugVolumes ? ALWAYS : outputLevel;
0920   if ( mot_vol->IsA() != TGeoVolumeAssembly::Class() )    {
0921     return nullptr;
0922   }
0923   Volume _v(mot_vol);
0924   if ( _v.testFlagBit(Volume::VETO_SIMU) )  {
0925     printout(lvl, "Geant4Converter", "++ AssemblyNode %s not converted [Veto'ed for simulation]",node->GetName());
0926     return nullptr;
0927   }
0928   Geant4GeometryInfo& info = data();
0929   Geant4AssemblyVolume* g4 = info.g4AssemblyVolumes[node];
0930   if ( g4 )  {
0931     printout(ALWAYS, "Geant4Converter", "+++ Assembly: **** : Re-using existing assembly: %s",node->GetName());
0932   }
0933   if ( !g4 )  {
0934     g4 = new Geant4AssemblyVolume();
0935     for(Int_t i=0; i < mot_vol->GetNdaughters(); ++i)   {
0936       TGeoNode*     dau     = mot_vol->GetNode(i);
0937       TGeoVolume*   dau_vol = dau->GetVolume();
0938       TGeoMatrix*   tr      = dau->GetMatrix();
0939       G4Transform3D transform;
0940 
0941       g4Transform(tr, transform);
0942       if ( is_left_handed(tr) )   {
0943         G4Scale3D     scale;
0944         G4Rotate3D    rot;
0945         G4Translate3D trans;
0946         transform.getDecomposition(scale, rot, trans);
0947         printout(debugReflections ? ALWAYS : lvl, "Geant4Converter",
0948                  "++ Placing reflected ASSEMBLY. dau:%s to mother %s "
0949                  "Tr:x=%8.1f y=%8.1f z=%8.1f   Scale:x=%4.2f y=%4.2f z=%4.2f",
0950                  dau_vol->GetName(), mot_vol->GetName(),
0951                  transform.dx(), transform.dy(), transform.dz(),
0952                  scale.xx(), scale.yy(), scale.zz());
0953       }
0954 
0955       if ( dau_vol->IsA() == TGeoVolumeAssembly::Class() )  {
0956         Geant4GeometryMaps::AssemblyMap::iterator ia = info.g4AssemblyVolumes.find(dau);
0957         if ( ia == info.g4AssemblyVolumes.end() )  {
0958           printout(FATAL, "Geant4Converter", "+++ Invalid child assembly at %s : %d  parent: %s child:%s",
0959                    __FILE__, __LINE__, name.c_str(), dau->GetName());
0960           delete g4;
0961           return nullptr;
0962         }
0963         g4->placeAssembly(dau, (*ia).second, transform);
0964         printout(lvl, "Geant4Converter", "+++ Assembly: AddPlacedAssembly %p: dau:%s "
0965                  "to mother %s Tr:x=%8.3f y=%8.3f z=%8.3f",
0966                  (void*)dau_vol, dau_vol->GetName(), mot_vol->GetName(),
0967                  transform.dx(), transform.dy(), transform.dz());
0968       }
0969       else   {
0970         Geant4GeometryMaps::VolumeMap::iterator iv = info.g4Volumes.find(dau_vol);
0971         if ( iv == info.g4Volumes.end() )  {
0972           printout(FATAL,"Geant4Converter", "+++ Invalid child volume at %s : %d  parent: %s child:%s",
0973                    __FILE__, __LINE__, name.c_str(), dau->GetName());
0974           except("Geant4Converter", "+++ Invalid child volume at %s : %d  parent: %s child:%s",
0975                  __FILE__, __LINE__, name.c_str(), dau->GetName());
0976         }
0977         g4->placeVolume(dau,(*iv).second, transform);
0978         printout(lvl, "Geant4Converter", "+++ Assembly: AddPlacedVolume %p: dau:%s "
0979                  "to mother %s Tr:x=%8.3f y=%8.3f z=%8.3f",
0980                  (void*)dau_vol, dau_vol->GetName(), mot_vol->GetName(),
0981                  transform.dx(), transform.dy(), transform.dz());
0982       }
0983     }
0984     info.g4AssemblyVolumes[node] = g4;
0985   }
0986   return g4;
0987 }
0988 
0989 /// Dump volume placement in GDML format to output stream
0990 void* Geant4Converter::handlePlacement(const std::string& name, const TGeoNode* node) const {
0991   Geant4GeometryInfo& info = this->data();
0992   PrintLevel lvl = this->debugPlacements ? ALWAYS : this->outputLevel;
0993   Geant4GeometryMaps::PlacementMap::const_iterator g4it = info.g4Placements.find(node);
0994   G4VPhysicalVolume* g4 = (g4it == info.g4Placements.end()) ? 0 : (*g4it).second;
0995   TGeoVolume* vol = node->GetVolume();
0996   Volume _v(vol);
0997 
0998   if ( _v.testFlagBit(Volume::VETO_SIMU) )  {
0999     printout(lvl, "Geant4Converter", "++ Placement %s not converted [Veto'ed for simulation]",node->GetName());
1000     return nullptr;
1001   }
1002   //g4 = nullptr;
1003   if ( !g4 ) {
1004     TGeoVolume* mot_vol = node->GetMotherVolume();
1005     TGeoMatrix* tr = node->GetMatrix();
1006     if ( !tr ) {
1007       except("Geant4Converter",
1008              "++ Attempt to handle placement without transformation:%p %s of type %s vol:%p",
1009              node, node->GetName(), node->IsA()->GetName(), vol);
1010     }
1011     else if (nullptr == vol) {
1012       except("Geant4Converter", "++ Unknown G4 volume:%p %s of type %s ptr:%p",
1013              node, node->GetName(), node->IsA()->GetName(), vol);
1014     }
1015     else {
1016       int  copy               = node->GetNumber();
1017       bool node_is_reflected  = is_left_handed(tr);
1018       bool node_is_assembly   = vol->IsA() == TGeoVolumeAssembly::Class();
1019       bool mother_is_assembly = mot_vol ? mot_vol->IsA() == TGeoVolumeAssembly::Class() : false;
1020 
1021       if ( mother_is_assembly )   {
1022         //
1023         // Mother is an assembly:
1024         // Nothing to do here, because:
1025         // -- placed volumes were already added before in "handleAssembly"
1026         // -- imprint cannot be made, because this requires a logical volume as a mother
1027         //
1028         printout(lvl, "Geant4Converter", "+++ Assembly: **** : daughter %s to mother %s",
1029                  vol->GetName(), mot_vol ? mot_vol->GetName() : "????");
1030         return nullptr;
1031       }
1032 
1033       G4LogicalVolume* g4mot = nullptr;
1034       auto             volIt = info.g4Volumes.find(mot_vol);
1035       if ( volIt != info.g4Volumes.end() )  {
1036         // Rational: Conversion/scanning was steered from top->down.
1037         //
1038         // Hence all mothers have proper relationship also to re-aligned daughters.
1039         // However, only the re-aligned daughters have the proper mother-relationship.
1040         // --> If a mother was found, the converted mother is the real one.
1041         g4mot = (*volIt).second;
1042       }
1043       else if ( node != info.manager->GetTopNode() )  {
1044         /// The mother of daughter was not converted.
1045         /// This may be a non-re-aligned daughter, where the daughter->mother relationship
1046         /// was not properly updated on re-alignment.
1047         ///
1048         TGeoNode *n1;
1049         TGeoIterator iter(info.manager->GetTopVolume());
1050         iter.SetType(1);
1051         printout(ALWAYS, "Geant4Converter", "+++ (SHOULD NOT ENTER HERE) Assembly: no G4 mother: %s org mot: %p",
1052                  node->GetName(), mot_vol);
1053         while ( (n1=iter.Next()) )  {
1054           if ( n1 == node )  {
1055             TGeoNode*   nmot = iter.GetNode(iter.GetLevel()-1);
1056             TGeoVolume* mmot = nmot->GetVolume();
1057             volIt = info.g4Volumes.find(mmot);
1058             if ( volIt != info.g4Volumes.end() )  {
1059               TString path;
1060               iter.GetPath(path);
1061               g4mot = (*volIt).second;
1062               printout(ALWAYS, "Geant4Converter", "+++ Assembly: Realigned mother: %s org mot: %p aligned: %p",
1063                        path.Data(), mot_vol, mmot);
1064               break;
1065             }
1066           }
1067         }
1068       }
1069 
1070       G4Scale3D        scale;
1071       G4Rotate3D       rotate;
1072       G4Translate3D    trans;
1073       G4Transform3D    transform;
1074       g4Transform(tr, transform);
1075       transform.getDecomposition(scale, rotate, trans);
1076       if ( node_is_assembly )   {
1077         //
1078         // Node is an assembly:
1079         // Imprint the assembly. The mother MUST already be transformed.
1080         //
1081         printout(lvl, "Geant4Converter", "++ Assembly: makeImprint: dau:%-12s %s in mother %-12s "
1082                  "Tr:x=%8.1f y=%8.1f z=%8.1f   Scale:x=%4.2f y=%4.2f z=%4.2f",
1083                  node->GetName(), node_is_reflected ? "(REFLECTED)" : "",
1084                  mot_vol ? mot_vol->GetName() : "<unknown>",
1085                  transform.dx(), transform.dy(), transform.dz(),
1086                  scale.xx(), scale.yy(), scale.zz());
1087         Geant4AssemblyVolume* ass = info.g4AssemblyVolumes[node];
1088         Geant4AssemblyVolume::Chain chain;
1089         chain.emplace_back(node);
1090         if ( !ass )  {
1091           except("Geant4Converter",
1092                  "+++ Assembly: %s mother: %s Geant4AssemblyVolume not present!",
1093                  node->GetName(), mot_vol ? mot_vol->GetName() : "<unknown>");
1094         }
1095         ass->imprint(*this, node, chain, ass, g4mot, transform, checkOverlaps);
1096         return nullptr;
1097       }
1098       else if ( node != info.manager->GetTopNode() && nullptr == g4mot )  {
1099         //throw std::logic_error("Geant4Converter: Invalid mother volume found!");
1100       }
1101       PlacedVolume pv(node);
1102       const auto*  pv_data = pv.data();
1103       G4LogicalVolume* g4vol = info.g4Volumes[vol];
1104       //G4LogicalVolume* g4mot = info.g4Volumes[mot_vol];
1105       G4PhysicalVolumesPair pvPlaced  { nullptr, nullptr };
1106 
1107       if ( pv_data && pv_data->params && (pv_data->params->flags&Volume::REPLICATED) )   {
1108         EAxis  axis = kUndefined;
1109         double width = 0e0, offset = 0e0;
1110         auto flags = pv_data->params->flags;
1111         auto count = pv_data->params->trafo1D.second;
1112         auto start = pv_data->params->start.Translation().Vect();
1113         auto delta = pv_data->params->trafo1D.first.Translation().Vect();
1114 
1115         if ( flags&Volume::X_axis )
1116         { axis = kXAxis; width = delta.X(); offset = start.X(); }
1117         else if ( flags&Volume::Y_axis )
1118         { axis = kYAxis; width = delta.Y(); offset = start.Y(); }
1119         else if ( flags&Volume::Z_axis )
1120         { axis = kZAxis; width = delta.Z(); offset = start.Z(); }
1121         else
1122           except("Geant4Converter",
1123                  "++ Replication around unknown axis is not implemented. flags: %16X", flags);
1124         printout(INFO,"Geant4Converter","++ Replicate: Axis: %ld Count: %ld offset: %f width: %f",
1125                  axis, count, offset, width);
1126         auto* g4pv = new G4PVReplica(name,      // its name
1127                                      g4vol,     // its logical volume
1128                                      g4mot,     // its mother (logical) volume
1129                                      axis,      // its replication axis
1130                                      count,     // Number of replicas
1131                                      width,     // Distance between 2 replicas
1132                                      offset);   // Placement offset in axis direction
1133         pvPlaced = { g4pv, nullptr };
1134 #if 0
1135         pvPlaced =
1136           G4ReflectionFactory::Instance()->Replicate(name,      // its name
1137                                                      g4vol,     // its logical volume
1138                                                      g4mot,     // its mother (logical) volume
1139                                                      axis,      // its replication axis
1140                                                      count,     // Number of replicas
1141                                                      width,     // Distance between 2 replicas
1142                                                      offset);   // Placement offset in axis direction
1143         /// Update replica list to avoid additional conversions...
1144         auto* g4pv = pvPlaced.second ? pvPlaced.second : pvPlaced.first;
1145 #endif
1146         for( auto& handle : pv_data->params->placements )
1147           info.g4Placements[handle.ptr()] = g4pv;
1148       }
1149       else if ( pv_data && pv_data->params )   {
1150         auto*  g4par = new Geant4PlacementParameterisation(pv);
1151         auto*  g4pv  = new G4PVParameterised(name,              // its name
1152                                              g4vol,             // its logical volume
1153                                              g4mot,             // its mother (logical) volume
1154                                              g4par->axis(),     // its replication axis
1155                                              g4par->count(),    // Number of replicas
1156                                              g4par);            // G4 parametrization
1157         pvPlaced = { g4pv, nullptr };
1158         /// Update replica list to avoid additional conversions...
1159         for( auto& handle : pv_data->params->placements )
1160           info.g4Placements[handle.ptr()] = g4pv;
1161       }
1162       else    {
1163         pvPlaced =
1164           G4ReflectionFactory::Instance()->Place(transform,     // no rotation
1165                                                  name,          // its name
1166                                                  g4vol,         // its logical volume
1167                                                  g4mot,         // its mother (logical) volume
1168                                                  false,         // no boolean operations
1169                                                  copy,          // its copy number
1170                                                  checkOverlaps);
1171       }
1172       printout(debugReflections||debugPlacements ? ALWAYS : lvl, "Geant4Converter",
1173                "++ Place %svolume %-12s in mother %-12s "
1174                "Tr:x=%8.1f y=%8.1f z=%8.1f   Scale:x=%4.2f y=%4.2f z=%4.2f",
1175                node_is_reflected ? "REFLECTED " : "", _v.name(),
1176                mot_vol ? mot_vol->GetName() : "<unknown>",
1177                transform.dx(), transform.dy(), transform.dz(),
1178                scale.xx(), scale.yy(), scale.zz());
1179       // First 2 cases can be combined.
1180       // Leave them separated for debugging G4ReflectionFactory for now...
1181       if ( node_is_reflected  && !pvPlaced.second )
1182         return info.g4Placements[node] = pvPlaced.first;
1183       else if ( !node_is_reflected && !pvPlaced.second )
1184         return info.g4Placements[node] = pvPlaced.first;
1185       // Now deal with valid pvPlaced.second ...
1186       if ( node_is_reflected )
1187         return info.g4Placements[node] = pvPlaced.first;
1188       else if ( !node_is_reflected )
1189         return info.g4Placements[node] = pvPlaced.first;
1190       g4 = pvPlaced.second ? pvPlaced.second : pvPlaced.first;
1191     }
1192     info.g4Placements[node] = g4;
1193     printout(ERROR, "Geant4Converter", "++ DEAD code. Should not end up here!");
1194   }
1195   return g4;
1196 }
1197 
1198 /// Convert the geometry type region into the corresponding Geant4 object(s).
1199 void* Geant4Converter::handleRegion(Region region, const std::set<const TGeoVolume*>& /* volumes */) const {
1200   G4Region* g4 = data().g4Regions[region];
1201   if ( !g4 ) {
1202     PrintLevel lvl = debugRegions ? ALWAYS : outputLevel;
1203     Region r = region;
1204     g4 = new G4Region(region.name());
1205 
1206     // create region info with storeSecondaries flag
1207     if( not r.wasThresholdSet() and r.storeSecondaries() ) {
1208       throw std::runtime_error("G4Region: StoreSecondaries is True, but no explicit threshold set:");
1209     }
1210     printout(lvl, "Geant4Converter", "++ Setting up region: %s", r.name());
1211     G4UserRegionInformation* info = new G4UserRegionInformation();
1212     info->region = r;
1213     info->threshold = r.threshold()*CLHEP::MeV/units::MeV;
1214     info->storeSecondaries = r.storeSecondaries();
1215     g4->SetUserInformation(info);
1216 
1217     printout(lvl, "Geant4Converter", "++ Converted region settings of:%s.", r.name());
1218     std::vector < std::string > &limits = r.limits();
1219     G4ProductionCuts* cuts = 0;
1220     // set production cut
1221     if( not r.useDefaultCut() ) {
1222       cuts = new G4ProductionCuts();
1223       cuts->SetProductionCut(r.cut()*CLHEP::mm/units::mm);
1224       printout(lvl, "Geant4Converter", "++ %s: Using default cut: %f [mm]",
1225                r.name(), r.cut()*CLHEP::mm/units::mm);
1226     }
1227     for( const auto& nam : limits )  {
1228       LimitSet ls = m_detDesc.limitSet(nam);
1229       if ( ls.isValid() ) {
1230         const LimitSet::Set& cts = ls.cuts();
1231         for (const auto& c : cts )   {
1232           int pid = 0;
1233           if ( c.particles == "*" ) pid = -1;
1234           else if ( c.particles == "e-"     ) pid = idxG4ElectronCut;
1235           else if ( c.particles == "e+"     ) pid = idxG4PositronCut;
1236           else if ( c.particles == "e[+-]"  ) pid = -idxG4PositronCut-idxG4ElectronCut;
1237           else if ( c.particles == "e[-+]"  ) pid = -idxG4PositronCut-idxG4ElectronCut;
1238           else if ( c.particles == "gamma"  ) pid = idxG4GammaCut;
1239           else if ( c.particles == "proton" ) pid = idxG4ProtonCut;
1240           else throw std::runtime_error("G4Region: Invalid production cut particle-type:" + c.particles);
1241           if ( !cuts ) cuts = new G4ProductionCuts();
1242           if ( pid == -(idxG4PositronCut+idxG4ElectronCut) )  {
1243             cuts->SetProductionCut(c.value*CLHEP::mm/units::mm, idxG4PositronCut);
1244             cuts->SetProductionCut(c.value*CLHEP::mm/units::mm, idxG4ElectronCut);
1245           }
1246           else  {
1247             cuts->SetProductionCut(c.value*CLHEP::mm/units::mm, pid);
1248           }
1249           printout(lvl, "Geant4Converter", "++ %s: Set cut  [%s/%d] = %f [mm]",
1250                    r.name(), c.particles.c_str(), pid, c.value*CLHEP::mm/units::mm);
1251         }
1252         bool found = false;
1253         const auto& lm = data().g4Limits;
1254         for (const auto& j : lm )   {
1255           if (nam == j.first->GetName()) {
1256             g4->SetUserLimits(j.second);
1257             printout(lvl, "Geant4Converter", "++ %s: Set limits %s to region type %s",
1258                      r.name(), nam.c_str(), j.second->GetType().c_str());
1259             found = true;
1260             break;
1261           }
1262         }
1263         if ( found )   {
1264           continue;
1265         }
1266       }
1267       except("Geant4Converter", "++ G4Region: Failed to resolve limitset: " + nam);
1268     }
1269     /// Assign cuts to region if they were created
1270     if ( cuts ) g4->SetProductionCuts(cuts);
1271     data().g4Regions[region] = g4;
1272   }
1273   return g4;
1274 }
1275 
1276 /// Convert the geometry type LimitSet into the corresponding Geant4 object(s).
1277 void* Geant4Converter::handleLimitSet(LimitSet limitset, const std::set<const TGeoVolume*>& /* volumes */) const {
1278   G4UserLimits* g4 = data().g4Limits[limitset];
1279   if ( !g4 ) {
1280     PrintLevel lvl = debugLimits || debugRegions ? ALWAYS : outputLevel;
1281     struct LimitPrint  {
1282       const LimitSet& ls;
1283       LimitPrint(const LimitSet& lset) : ls(lset) {}
1284       const LimitPrint& operator()(const std::string& pref, const Geant4UserLimits::Handler& h)  const {
1285         if ( !h.particleLimits.empty() )  {
1286           printout(ALWAYS,"Geant4Converter",
1287                    "+++ LimitSet: Explicit Limit %s.%s applied for particles:",ls.name(), pref.c_str());
1288           for(const auto& p : h.particleLimits)
1289             printout(ALWAYS,"Geant4Converter","+++ LimitSet:    Particle type: %-18s PDG: %-6d : %f",
1290                      p.first->GetParticleName().c_str(), p.first->GetPDGEncoding(), p.second);
1291         }
1292         else if ( h.defaultValue > std::numeric_limits<double>::epsilon() )  {
1293           printout(ALWAYS,"Geant4Converter",
1294                    "+++ LimitSet: Implicit Limit %s.%s for wildcard particles: %f",
1295                    ls.name(), pref.c_str(), float(h.defaultValue));
1296         }
1297         return *this;
1298       }
1299     };
1300     Geant4UserLimits* limits = new Geant4UserLimits(limitset);
1301     g4 = limits;
1302     printout(lvl, "Geant4Converter",
1303              "++ Successfully converted LimitSet: %s [%ld cuts, %ld limits]",
1304              limitset.name(), limitset.cuts().size(), limitset.limits().size());
1305     if ( debugRegions || debugLimits )    {
1306       LimitPrint print(limitset);
1307       print("maxTime",    limits->maxTime)
1308         ("minEKine",      limits->minEKine)
1309         ("minRange",      limits->minRange)
1310         ("maxStepLength", limits->maxStepLength)
1311         ("maxTrackLength",limits->maxTrackLength);
1312     }
1313     data().g4Limits[limitset] = g4;
1314   }
1315   return g4;
1316 }
1317 
1318 /// Convert the geometry visualisation attributes to the corresponding Geant4 object(s).
1319 void* Geant4Converter::handleVis(const std::string& /* name */, VisAttr attr) const {
1320   Geant4GeometryInfo& info = data();
1321   G4VisAttributes*    g4   = info.g4Vis[attr];
1322   if ( !g4 ) {
1323     float red = 0, green = 0, blue = 0;
1324     int   style = attr.lineStyle();
1325     attr.rgb(red, green, blue);
1326     g4 = new G4VisAttributes(attr.visible(), G4Colour(red, green, blue, attr.alpha()));
1327     //g4->SetLineWidth(attr->GetLineWidth());
1328     g4->SetDaughtersInvisible(!attr.showDaughters());
1329     if ( style == VisAttr::SOLID ) {
1330       g4->SetLineStyle(G4VisAttributes::unbroken);
1331       g4->SetForceWireframe(false);
1332       g4->SetForceSolid(true);
1333     }
1334     else if ( style == VisAttr::WIREFRAME || style == VisAttr::DASHED ) {
1335       g4->SetLineStyle(G4VisAttributes::dashed);
1336       g4->SetForceSolid(false);
1337       g4->SetForceWireframe(true);
1338     }
1339     info.g4Vis[attr] = g4;
1340   }
1341   return g4;
1342 }
1343 
1344 /// Handle the geant 4 specific properties
1345 void Geant4Converter::handleProperties(Detector::Properties& prp) const {
1346   std::map < std::string, std::string > processors;
1347   static int s_idd = 9999999;
1348   for( const auto& [nam, vals] : prp ) {
1349     if ( nam.substr(0, 6) == "geant4" ) {
1350       auto id_it = vals.find("id");
1351       std::string id = (id_it == vals.end()) ? _toString(++s_idd,"%d") : (*id_it).second;
1352       processors.emplace(id, nam);
1353     }
1354   }
1355   for( const auto& p : processors ) {
1356     const GeoHandler* hdlr = this;
1357     const Detector::PropertyValues& vals = prp[p.second];
1358     auto iter = vals.find("type");
1359     if ( iter != vals.end() )  {
1360       std::string type = iter->second;
1361       std::string tag  = type + "_Geant4_action";
1362       Detector* det = const_cast<Detector*>(&m_detDesc);
1363       long      res = PluginService::Create<long>(tag, det, hdlr, &vals);
1364       if ( 0 == res ) {
1365         throw std::runtime_error("Failed to locate plugin to interprete files of type"
1366                                  " \"" + tag + "\" - no factory:" + type);
1367       }
1368       res = *(long*)res;
1369       if ( res != 1 ) {
1370         throw std::runtime_error("Failed to invoke the plugin " + tag + " of type " + type);
1371       }
1372       printout(outputLevel, "Geant4Converter", "+++++ Executed Successfully Geant4 setup module *%s*.", type.c_str());
1373       continue;
1374     }
1375     printout(outputLevel, "Geant4Converter", "+++++ FAILED to execute Geant4 setup module *%s*.", p.second.c_str());    
1376   }
1377 }
1378 
1379 /// Convert the geometry type material into the corresponding Geant4 object(s).
1380 void* Geant4Converter::handleMaterialProperties(TObject* mtx) const    {
1381   Geant4GeometryInfo& info   = data();
1382   TGDMLMatrix*        matrix = (TGDMLMatrix*)mtx;
1383   const char*         cptr   = ::strstr(matrix->GetName(), GEANT4_TAG_IGNORE);
1384   Geant4GeometryInfo::PropertyVector* g4 = info.g4OpticalProperties[matrix];
1385 
1386   if ( nullptr != cptr )   {  // Check if the property should not be passed to Geant4
1387     printout(INFO,"Geant4MaterialProperties","++ Ignore property %s [%s].",
1388              matrix->GetName(), matrix->GetTitle());             
1389     return nullptr;
1390   }
1391   cptr = ::strstr(matrix->GetTitle(), GEANT4_TAG_IGNORE);
1392   if ( nullptr != cptr )   {  // Check if the property should not be passed to Geant4
1393     printout(INFO,"Geant4MaterialProperties","++ Ignore property %s [%s].",
1394              matrix->GetName(), matrix->GetTitle());
1395     return nullptr;
1396   }
1397   
1398   if ( !g4 )  {
1399     PrintLevel lvl = debugMaterials ? ALWAYS : outputLevel;
1400     g4 = new Geant4GeometryInfo::PropertyVector();
1401     std::size_t rows = matrix->GetRows();
1402     g4->name    = matrix->GetName();
1403     g4->title   = matrix->GetTitle();
1404     g4->bins.reserve(rows);
1405     g4->values.reserve(rows);
1406     for( std::size_t i=0; i<rows; ++i )   {
1407       g4->bins.emplace_back(matrix->Get(i,0)  /*   *CLHEP::eV/units::eV   */);
1408       g4->values.emplace_back(matrix->Get(i,1));
1409     }
1410     printout(lvl, "Geant4Converter",
1411              "++ Successfully converted material property:%s : %s [%ld rows]",
1412              matrix->GetName(), matrix->GetTitle(), rows);
1413     info.g4OpticalProperties[matrix] = g4;
1414   }
1415   return g4;
1416 }
1417 
1418 static G4OpticalSurfaceFinish geant4_surface_finish(TGeoOpticalSurface::ESurfaceFinish f)   {
1419 #define TO_G4_FINISH(x)  case TGeoOpticalSurface::kF##x : return x;
1420   switch(f)   {
1421     TO_G4_FINISH(polished);              // smooth perfectly polished surface
1422     TO_G4_FINISH(polishedfrontpainted);  // smooth top-layer (front) paint
1423     TO_G4_FINISH(polishedbackpainted);   // same is 'polished' but with a back-paint
1424  
1425     TO_G4_FINISH(ground);                // rough surface
1426     TO_G4_FINISH(groundfrontpainted);    // rough top-layer (front) paint
1427     TO_G4_FINISH(groundbackpainted);     // same as 'ground' but with a back-paint
1428 
1429     TO_G4_FINISH(polishedlumirrorair);   // mechanically polished surface, with lumirror
1430     TO_G4_FINISH(polishedlumirrorglue);  // mechanically polished surface, with lumirror & meltmount
1431     TO_G4_FINISH(polishedair);           // mechanically polished surface
1432     TO_G4_FINISH(polishedteflonair);     // mechanically polished surface, with teflon
1433     TO_G4_FINISH(polishedtioair);        // mechanically polished surface, with tio paint
1434     TO_G4_FINISH(polishedtyvekair);      // mechanically polished surface, with tyvek
1435     TO_G4_FINISH(polishedvm2000air);     // mechanically polished surface, with esr film
1436     TO_G4_FINISH(polishedvm2000glue);    // mechanically polished surface, with esr film & meltmount
1437 
1438     TO_G4_FINISH(etchedlumirrorair);     // chemically etched surface, with lumirror
1439     TO_G4_FINISH(etchedlumirrorglue);    // chemically etched surface, with lumirror & meltmount
1440     TO_G4_FINISH(etchedair);             // chemically etched surface
1441     TO_G4_FINISH(etchedteflonair);       // chemically etched surface, with teflon
1442     TO_G4_FINISH(etchedtioair);          // chemically etched surface, with tio paint
1443     TO_G4_FINISH(etchedtyvekair);        // chemically etched surface, with tyvek
1444     TO_G4_FINISH(etchedvm2000air);       // chemically etched surface, with esr film
1445     TO_G4_FINISH(etchedvm2000glue);      // chemically etched surface, with esr film & meltmount
1446 
1447     TO_G4_FINISH(groundlumirrorair);     // rough-cut surface, with lumirror
1448     TO_G4_FINISH(groundlumirrorglue);    // rough-cut surface, with lumirror & meltmount
1449     TO_G4_FINISH(groundair);             // rough-cut surface
1450     TO_G4_FINISH(groundteflonair);       // rough-cut surface, with teflon
1451     TO_G4_FINISH(groundtioair);          // rough-cut surface, with tio paint
1452     TO_G4_FINISH(groundtyvekair);        // rough-cut surface, with tyvek
1453     TO_G4_FINISH(groundvm2000air);       // rough-cut surface, with esr film
1454     TO_G4_FINISH(groundvm2000glue);      // rough-cut surface, with esr film & meltmount
1455 
1456     // for DAVIS model
1457     TO_G4_FINISH(Rough_LUT);             // rough surface
1458     TO_G4_FINISH(RoughTeflon_LUT);       // rough surface wrapped in Teflon tape
1459     TO_G4_FINISH(RoughESR_LUT);          // rough surface wrapped with ESR
1460     TO_G4_FINISH(RoughESRGrease_LUT);    // rough surface wrapped with ESR and coupled with opical grease
1461     TO_G4_FINISH(Polished_LUT);          // polished surface
1462     TO_G4_FINISH(PolishedTeflon_LUT);    // polished surface wrapped in Teflon tape
1463     TO_G4_FINISH(PolishedESR_LUT);       // polished surface wrapped with ESR
1464     TO_G4_FINISH(PolishedESRGrease_LUT); // polished surface wrapped with ESR and coupled with opical grease
1465     TO_G4_FINISH(Detector_LUT);          // polished surface with optical grease
1466   default:
1467     printout(ERROR,"Geant4Surfaces","++ Unknown finish style: %d [%s]. Assume polished!",
1468              int(f), TGeoOpticalSurface::FinishToString(f));
1469     return polished;
1470   }
1471 #undef TO_G4_FINISH
1472 }
1473 
1474 static G4SurfaceType geant4_surface_type(TGeoOpticalSurface::ESurfaceType t)   {
1475 #define TO_G4_TYPE(x)  case TGeoOpticalSurface::kT##x : return x;
1476   switch(t)   {
1477     TO_G4_TYPE(dielectric_metal);      // dielectric-metal interface
1478     TO_G4_TYPE(dielectric_dielectric); // dielectric-dielectric interface
1479     TO_G4_TYPE(dielectric_LUT);        // dielectric-Look-Up-Table interface
1480     TO_G4_TYPE(dielectric_LUTDAVIS);   // dielectric-Look-Up-Table DAVIS interface
1481     TO_G4_TYPE(dielectric_dichroic);   // dichroic filter interface
1482     TO_G4_TYPE(firsov);                // for Firsov Process
1483     TO_G4_TYPE(x_ray);                  // for x-ray mirror process
1484   default:
1485     printout(ERROR,"Geant4Surfaces","++ Unknown surface type: %d [%s]. Assume dielectric_metal!",
1486              int(t), TGeoOpticalSurface::TypeToString(t));
1487     return dielectric_metal;
1488   }
1489 #undef TO_G4_TYPE
1490 }
1491 
1492 static G4OpticalSurfaceModel geant4_surface_model(TGeoOpticalSurface::ESurfaceModel surfMod)   {
1493 #define TO_G4_MODEL(x)  case TGeoOpticalSurface::kM##x : return x;
1494   switch(surfMod)   {
1495     TO_G4_MODEL(glisur);   // original GEANT3 model
1496     TO_G4_MODEL(unified);  // UNIFIED model
1497     TO_G4_MODEL(LUT);      // Look-Up-Table model
1498     TO_G4_MODEL(DAVIS);    // DAVIS model
1499     TO_G4_MODEL(dichroic); // dichroic filter
1500   default:
1501     printout(ERROR,"Geant4Surfaces","++ Unknown surface model: %d [%s]. Assume glisur!",
1502              int(surfMod), TGeoOpticalSurface::ModelToString(surfMod));
1503     return glisur;
1504   }
1505 #undef TO_G4_MODEL
1506 }
1507 
1508 /// Convert the optical surface to Geant4
1509 void* Geant4Converter::handleOpticalSurface(TObject* surface) const    {
1510   TGeoOpticalSurface* optSurf    = (TGeoOpticalSurface*)surface;
1511   Geant4GeometryInfo& info = data();
1512   G4OpticalSurface*   g4   = info.g4OpticalSurfaces[optSurf];
1513   if ( !g4 ) {
1514     G4SurfaceType          type   = geant4_surface_type(optSurf->GetType());
1515     G4OpticalSurfaceModel  model  = geant4_surface_model(optSurf->GetModel());
1516     G4OpticalSurfaceFinish finish = geant4_surface_finish(optSurf->GetFinish());
1517     std::string name = make_NCName(optSurf->GetName());
1518     PrintLevel lvl = debugSurfaces ? ALWAYS : DEBUG;
1519     g4 = new G4OpticalSurface(name, model, finish, type, optSurf->GetValue());
1520     g4->SetSigmaAlpha(optSurf->GetSigmaAlpha());
1521     g4->SetPolish(optSurf->GetPolish());
1522 
1523     printout(lvl, "Geant4Converter",
1524              "++ Created OpticalSurface: %-18s type:%s model:%s finish:%s SigmaAlphs: %.3e Polish: %.3e",
1525              optSurf->GetName(),
1526              TGeoOpticalSurface::TypeToString(optSurf->GetType()),
1527              TGeoOpticalSurface::ModelToString(optSurf->GetModel()),
1528              TGeoOpticalSurface::FinishToString(optSurf->GetFinish()),
1529              optSurf->GetSigmaAlpha(), optSurf->GetPolish());
1530     ///
1531     /// Convert non-scalar properties from GDML tables
1532     G4MaterialPropertiesTable* tab = nullptr;
1533     TListIter itp(&optSurf->GetProperties());
1534     for(TObject* obj = itp.Next(); obj; obj = itp.Next())  {
1535       std::string exc_str;
1536       TNamed*      named  = (TNamed*)obj;
1537       TGDMLMatrix* matrix = info.manager->GetGDMLMatrix(named->GetTitle());
1538       const char*  cptr   = ::strstr(matrix->GetName(), GEANT4_TAG_IGNORE);
1539       if ( nullptr != cptr )  // Check if the property should not be passed to Geant4
1540         continue;
1541 
1542       if ( nullptr == tab )  {
1543         tab = new G4MaterialPropertiesTable();
1544         g4->SetMaterialPropertiesTable(tab);
1545       }
1546 
1547       Geant4GeometryInfo::PropertyVector* v =
1548         (Geant4GeometryInfo::PropertyVector*)handleMaterialProperties(matrix);
1549       if ( !v )  {  // Error!
1550         except("Geant4OpticalSurface","++ Failed to convert opt.surface %s. Property table %s is not defined!",
1551                optSurf->GetName(), named->GetTitle());
1552       }
1553       int idx = -1;
1554       try   {
1555         idx = tab->GetPropertyIndex(named->GetName());
1556       }
1557       catch(const std::exception& e)   {
1558         exc_str = e.what();
1559       }
1560       catch(...)   {
1561       }
1562       if ( idx < 0 )   {
1563         printout(ERROR, "Geant4Converter",
1564                  "++ UNKNOWN Geant4 Property: %-20s %s [IGNORED]",
1565                  exc_str.c_str(), named->GetName());
1566         continue;
1567       }
1568       // We need to convert the property from TGeo units to Geant4 units
1569       auto conv = g4PropertyConversion(idx);
1570       std::vector<double> bins(v->bins), vals(v->values);
1571       for(std::size_t i=0, count=v->bins.size(); i<count; ++i)
1572         bins[i] *= conv.first, vals[i] *= conv.second;
1573       G4MaterialPropertyVector* vec = new G4MaterialPropertyVector(&bins[0], &vals[0], bins.size());
1574       tab->AddProperty(named->GetName(), vec);
1575       
1576       printout(lvl, "Geant4Converter",
1577                "++       Property: %-20s [%ld x %ld] -->  %s",
1578                named->GetName(), matrix->GetRows(), matrix->GetCols(), named->GetTitle());
1579       for(std::size_t i=0, count=v->bins.size(); i<count; ++i)
1580         printout(lvl, named->GetName(),
1581                  "  Geant4: %8.3g [MeV]  TGeo: %8.3g [GeV] Conversion: %8.3g",
1582                  bins[i], v->bins[i], conv.first);
1583     }
1584     ///
1585     /// Convert scalar properties
1586 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,31,1)
1587     TListIter itc(&optSurf->GetConstProperties());
1588     for(TObject* obj = itc.Next(); obj; obj = itc.Next())  {
1589       std::string  exc_str;
1590       TNamed* named  = (TNamed*)obj;
1591       const char* cptr = ::strstr(named->GetName(), GEANT4_TAG_IGNORE);
1592       if ( nullptr != cptr )   {
1593         printout(INFO, name, "++ Ignore CONST property %s [%s].",
1594                  named->GetName(), named->GetTitle());
1595         continue;
1596       }
1597       cptr = ::strstr(named->GetTitle(), GEANT4_TAG_IGNORE);
1598       if ( nullptr != cptr )   {
1599         printout(INFO, name,"++ Ignore CONST property %s [%s].",
1600                  named->GetName(), named->GetTitle());
1601         continue;
1602       }
1603       Bool_t   err = kFALSE;
1604       Double_t value = info.manager->GetProperty(named->GetTitle(),&err);
1605       if ( err != kFALSE )   {
1606         except(name,
1607                "++ FAILED to create G4 material %s [Cannot convert const property: %s]",
1608                optSurf->GetName(), named->GetName());
1609       }
1610       if ( nullptr == tab )  {
1611         tab = new G4MaterialPropertiesTable();
1612         g4->SetMaterialPropertiesTable(tab);
1613       }
1614       int idx = -1;
1615       try   {
1616         idx = tab->GetConstPropertyIndex(named->GetName());
1617       }
1618       catch(const std::exception& e)   {
1619         exc_str = e.what();
1620       }
1621       catch(...)   {
1622       }
1623       if ( idx < 0 )   {
1624         printout(ERROR, name,
1625                  "++ UNKNOWN Geant4 CONST Property: %-20s %s [IGNORED]",
1626                  exc_str.c_str(), named->GetName());
1627         continue;
1628       }
1629       // We need to convert the property from TGeo units to Geant4 units
1630       double conv = g4ConstPropertyConversion(idx);
1631       printout(lvl, name, "++      CONST Property: %-20s %g * %g --> %g ",
1632                named->GetName(), value, conv, value * conv);
1633       tab->AddConstProperty(named->GetName(), value * conv);
1634     }
1635 #endif  // ROOT_VERSION >= 6.31.1
1636     info.g4OpticalSurfaces[optSurf] = g4;
1637   }
1638   return g4;
1639 }
1640 
1641 /// Convert the skin surface to Geant4
1642 void* Geant4Converter::handleSkinSurface(TObject* surface) const   {
1643   TGeoSkinSurface*    surf = (TGeoSkinSurface*)surface;
1644   Geant4GeometryInfo& info = data();
1645   G4LogicalSkinSurface* g4 = info.g4SkinSurfaces[surf];
1646   if ( !g4 ) {
1647     G4OpticalSurface* optSurf  = info.g4OpticalSurfaces[OpticalSurface(surf->GetSurface())];
1648     G4LogicalVolume*  v = info.g4Volumes[surf->GetVolume()];
1649     std::string name = make_NCName(surf->GetName());
1650     g4 = new G4LogicalSkinSurface(name, v, optSurf);
1651     printout(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter",
1652              "++ Created SkinSurface: %-18s  optical:%s",
1653              surf->GetName(), surf->GetSurface()->GetName());
1654     info.g4SkinSurfaces[surf] = g4;
1655   }
1656   return g4;
1657 }
1658 
1659 /// Convert the border surface to Geant4
1660 void* Geant4Converter::handleBorderSurface(TObject* surface) const   {
1661   TGeoBorderSurface*    surf = (TGeoBorderSurface*)surface;
1662   Geant4GeometryInfo&   info = data();
1663   G4LogicalBorderSurface* g4 = info.g4BorderSurfaces[surf];
1664   if ( !g4 ) {
1665     G4OpticalSurface*  optSurf = info.g4OpticalSurfaces[OpticalSurface(surf->GetSurface())];
1666     G4VPhysicalVolume* n1 = info.g4Placements[surf->GetNode1()];
1667     G4VPhysicalVolume* n2 = info.g4Placements[surf->GetNode2()];
1668     std::string name = make_NCName(surf->GetName());
1669     g4 = new G4LogicalBorderSurface(name, n1, n2, optSurf);
1670     printout(debugSurfaces ? ALWAYS : DEBUG, "Geant4Converter",
1671              "++ Created BorderSurface: %-18s  optical:%s",
1672              surf->GetName(), surf->GetSurface()->GetName());
1673     info.g4BorderSurfaces[surf] = g4;
1674   }
1675   return g4;
1676 }
1677 
1678 /// Convert the geometry type SensitiveDetector into the corresponding Geant4 object(s).
1679 void Geant4Converter::printSensitive(SensitiveDetector sens_det, const std::set<const TGeoVolume*>& /* volumes */) const {
1680   Geant4GeometryInfo&          info = data();
1681   std::set<const TGeoVolume*>& volset = info.sensitives[sens_det];
1682   SensitiveDetector            sd = sens_det;
1683   std::stringstream str;
1684 
1685   printout(INFO, "Geant4Converter", "++ SensitiveDetector: %-18s %-20s Hits:%-16s", sd.name(), ("[" + sd.type() + "]").c_str(),
1686            sd.hitsCollection().c_str());
1687   str << "                    | " << "Cutoff:" << std::setw(6) << std::left
1688       << sd.energyCutoff() << std::setw(5) << std::right << volset.size()
1689       << " volumes ";
1690   if (sd.region().isValid())
1691     str << " Region:" << std::setw(12) << std::left << sd.region().name();
1692   if (sd.limits().isValid())
1693     str << " Limits:" << std::setw(12) << std::left << sd.limits().name();
1694   str << ".";
1695   printout(INFO, "Geant4Converter", str.str().c_str());
1696 
1697   for (const auto i : volset )  {
1698     std::map<Volume, G4LogicalVolume*>::iterator v = info.g4Volumes.find(i);
1699     if ( v != info.g4Volumes.end() )   {
1700       G4LogicalVolume* vol = (*v).second;
1701       str.str("");
1702       str << "                                   | " << "Volume:" << std::setw(24) << std::left << vol->GetName() << " "
1703           << vol->GetNoDaughters() << " daughters.";
1704       printout(INFO, "Geant4Converter", str.str().c_str());
1705     }
1706   }
1707 }
1708 
1709 std::string printSolid(G4VSolid* sol) {
1710   std::stringstream str;
1711   if (typeid(*sol) == typeid(G4Box)) {
1712     const G4Box* b = (G4Box*) sol;
1713     str << "++ Box: x=" << b->GetXHalfLength() << " y=" << b->GetYHalfLength() << " z=" << b->GetZHalfLength();
1714   }
1715   else if (typeid(*sol) == typeid(G4Tubs)) {
1716     const G4Tubs* t = (const G4Tubs*) sol;
1717     str << " Tubs: Ri=" << t->GetInnerRadius() << " Ra=" << t->GetOuterRadius() << " z/2=" << t->GetZHalfLength() << " Phi="
1718         << t->GetStartPhiAngle() << "..." << t->GetDeltaPhiAngle();
1719   }
1720   return str.str();
1721 }
1722 
1723 /// Print G4 placement
1724 void* Geant4Converter::printPlacement(const std::string& name, const TGeoNode* node) const {
1725   Geant4GeometryInfo& info = data();
1726   G4VPhysicalVolume*  g4   = info.g4Placements[node];
1727   G4LogicalVolume*    vol  = info.g4Volumes[node->GetVolume()];
1728   G4LogicalVolume*    mot  = info.g4Volumes[node->GetMotherVolume()];
1729   G4VSolid*           sol  = vol->GetSolid();
1730   G4ThreeVector       tr   = g4->GetObjectTranslation();
1731   G4VSensitiveDetector* sd = vol->GetSensitiveDetector();
1732   if ( !sd )  {
1733     return g4;
1734   }
1735   std::stringstream str;
1736   str << "G4Cnv::placement: + " << name << " No:" << node->GetNumber() << " Vol:" << vol->GetName() << " Solid:"
1737       << sol->GetName();
1738   printout(outputLevel, "G4Placement", str.str().c_str());
1739   str.str("");
1740   str << "                  |" << " Loc: x=" << tr.x() << " y=" << tr.y() << " z=" << tr.z();
1741   printout(outputLevel, "G4Placement", str.str().c_str());
1742   printout(outputLevel, "G4Placement", printSolid(sol).c_str());
1743   str.str("");
1744   str << "                  |" << " Ndau:" << vol->GetNoDaughters()
1745       << " physvols." << " Mat:" << vol->GetMaterial()->GetName()
1746       << " Mother:" << (char*) (mot ? mot->GetName().c_str() : "---");
1747   printout(outputLevel, "G4Placement", str.str().c_str());
1748   str.str("");
1749   str << "                  |" << " SD:" << sd->GetName();
1750   printout(outputLevel, "G4Placement", str.str().c_str());
1751   return g4;
1752 }
1753 
1754 /// Create geometry conversion
1755 Geant4Converter& Geant4Converter::create(DetElement top) {
1756   typedef std::map<const TGeoNode*, std::vector<TGeoNode*> > _DAU;
1757   TTimeStamp start;
1758   _DAU daughters;
1759   Geant4GeometryInfo& geo = this->init();
1760   World wrld = top.world();
1761 
1762   m_data->clear();
1763   m_set_data->clear();
1764   m_daughters = &daughters;
1765   geo.manager = &wrld.detectorDescription().manager();
1766   this->collect(top, geo);
1767   this->checkOverlaps = false;
1768   // We do not have to handle defines etc.
1769   // All positions and the like are not really named.
1770   // Hence, start creating the G4 objects for materials, solids and log volumes.
1771   handleArray(this, geo.manager->GetListOfGDMLMatrices(), &Geant4Converter::handleMaterialProperties);
1772   handleArray(this, geo.manager->GetListOfOpticalSurfaces(), &Geant4Converter::handleOpticalSurface);
1773   
1774   handle(this,     geo.volumes, &Geant4Converter::collectVolume);
1775   handle(this,     geo.solids,  &Geant4Converter::handleSolid);
1776   printout(outputLevel, "Geant4Converter", "++ Handled %ld solids.", geo.solids.size());
1777   handleRefs(this, geo.vis,     &Geant4Converter::handleVis);
1778   printout(outputLevel, "Geant4Converter", "++ Handled %ld visualization attributes.", geo.vis.size());
1779   handleMap(this,  geo.limits,  &Geant4Converter::handleLimitSet);
1780   printout(outputLevel, "Geant4Converter", "++ Handled %ld limit sets.", geo.limits.size());
1781   handleMap(this,  geo.regions, &Geant4Converter::handleRegion);
1782   printout(outputLevel, "Geant4Converter", "++ Handled %ld regions.", geo.regions.size());
1783   handle(this,     geo.volumes, &Geant4Converter::handleVolume);
1784   printout(outputLevel, "Geant4Converter", "++ Handled %ld volumes.", geo.volumes.size());
1785   handleRMap(this, *m_data,     &Geant4Converter::handleAssembly);
1786   // Now place all this stuff appropriately
1787   //handleRMap(this, *m_data,     &Geant4Converter::handlePlacement);
1788   std::map<int, std::vector<const TGeoNode*> >::const_reverse_iterator i = m_data->rbegin();
1789   for ( ; i != m_data->rend(); ++i )  {
1790     for ( const TGeoNode* node : i->second )  {
1791       this->handlePlacement(node->GetName(), node);
1792     }
1793   }
1794   /// Handle concrete surfaces
1795   handleArray(this, geo.manager->GetListOfSkinSurfaces(),   &Geant4Converter::handleSkinSurface);
1796   handleArray(this, geo.manager->GetListOfBorderSurfaces(), &Geant4Converter::handleBorderSurface);
1797   //==================== Fields
1798   handleProperties(m_detDesc.properties());
1799   if ( printSensitives )  {
1800     handleMap(this, geo.sensitives, &Geant4Converter::printSensitive);
1801   }
1802   if ( printPlacements )  {
1803     handleRMap(this, *m_data, &Geant4Converter::printPlacement);
1804   }
1805 
1806   m_daughters = nullptr;
1807   geo.setWorld(top.placement().ptr());
1808   geo.valid = true;
1809   TTimeStamp stop;
1810   printout(INFO, "Geant4Converter",
1811            "+++  Successfully converted geometry to Geant4. [%7.3f seconds]",
1812            stop.AsDouble()-start.AsDouble() );
1813   return *this;
1814 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!