DDG4/src/Geant4Converter.cpp

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