Geant4Converter.cpp

//==========================================================================
//  AIDA Detector description implementation 
//--------------------------------------------------------------------------
// Copyright (C) Organisation europeenne pour la Recherche nucleaire (CERN)
// All rights reserved.
//
// For the licensing terms see $DD4hepINSTALL/LICENSE.
// For the list of contributors see $DD4hepINSTALL/doc/CREDITS.
//
// Author     : M.Frank
//
//==========================================================================

// Framework include files
#include <DD4hep/Shapes.h>
#include <DD4hep/Volumes.h>
#include <DD4hep/Plugins.h>
#include <DD4hep/Printout.h>
#include <DD4hep/Detector.h>
#include <DD4hep/DD4hepUnits.h>
#include <DD4hep/PropertyTable.h>
#include <DD4hep/DetectorTools.h>
#include <DD4hep/detail/ShapesInterna.h>
#include <DD4hep/detail/ObjectsInterna.h>
#include <DD4hep/detail/DetectorInterna.h>

#include <DDG4/Geant4Field.h>
#include <DDG4/Geant4Helpers.h>
#include <DDG4/Geant4Converter.h>
#include <DDG4/Geant4UserLimits.h>
#include <DDG4/Geant4AssemblyVolume.h>
#include <DDG4/Geant4PlacementParameterisation.h>
#include "Geant4ShapeConverter.h"

// ROOT includes
#include <TClass.h>
#include <TGeoBoolNode.h>

// Geant4 include files
#include <G4Version.hh>
#include <G4VisAttributes.hh>
#include <G4PVParameterised.hh>
#include <G4ProductionCuts.hh>
#include <G4VUserRegionInformation.hh>

#include <G4Box.hh>
#include <G4Tubs.hh>
#include <G4Ellipsoid.hh>
#include <G4UnionSolid.hh>
#include <G4ReflectedSolid.hh>
#include <G4SubtractionSolid.hh>
#include <G4IntersectionSolid.hh>
#include <G4VSensitiveDetector.hh>

#include <G4Region.hh>
#include <G4Element.hh>
#include <G4Isotope.hh>
#include <G4Material.hh>
#include <G4UserLimits.hh>
#include <G4RegionStore.hh>
#include <G4FieldManager.hh>
#include <G4LogicalVolume.hh>
#include <G4OpticalSurface.hh>
#include <G4ReflectionFactory.hh>
#include <G4LogicalSkinSurface.hh>
#include <G4ElectroMagneticField.hh>
#include <G4LogicalBorderSurface.hh>
#include <G4MaterialPropertiesTable.hh>
#if G4VERSION_NUMBER >= 1040
#include <G4MaterialPropertiesIndex.hh>
#endif
#include <G4ScaledSolid.hh>
#include <CLHEP/Units/SystemOfUnits.h>

// C/C++ include files
#include <iostream>
#include <iomanip>
#include <sstream>
#include <limits>

namespace units = dd4hep;
using namespace dd4hep::sim;
using namespace dd4hep;

namespace {

  static constexpr const double CM_2_MM = (CLHEP::centimeter/dd4hep::centimeter);
  static constexpr const char* GEANT4_TAG_IGNORE = "Geant4-ignore";
  static constexpr const char* GEANT4_TAG_PLUGIN = "Geant4-plugin";
  static constexpr const char* GEANT4_TAG_BIRKSCONSTANT    = "BirksConstant";
  static constexpr const char* GEANT4_TAG_MEE              = "MeanExcitationEnergy";
  static constexpr const char* GEANT4_TAG_ENE_PER_ION_PAIR = "MeanEnergyPerIonPair";

  static std::string indent = "";

  template <typename O, typename C, typename F> void handleRefs(const O* o, const C& c, F pmf) {
    for (typename C::const_iterator i = c.begin(); i != c.end(); ++i) {
      //(o->*pmf)((*i)->GetName(), *i);
      (o->*pmf)("", *i);
    }
  }

  template <typename O, typename C, typename F> void handle(const O* o, const C& c, F pmf) {
    for (typename C::const_iterator i = c.begin(); i != c.end(); ++i) {
      (o->*pmf)((*i)->GetName(), *i);
    }
  }

  template <typename O, typename F> void handleArray(const O* o, const TObjArray* c, F pmf) {
    TObjArrayIter arr(c);
    for(TObject* i = arr.Next(); i; i=arr.Next())
      (o->*pmf)(i);
  }

  template <typename O, typename C, typename F> void handleMap(const O* o, const C& c, F pmf) {
    for (typename C::const_iterator i = c.begin(); i != c.end(); ++i)
      (o->*pmf)((*i).first, (*i).second);
  }

  template <typename O, typename C, typename F> void handleRMap(const O* o, const C& c, F pmf) {
    for (typename C::const_reverse_iterator i = c.rbegin(); i != c.rend(); ++i)  {
      //cout << "Handle RMAP [ " << (*i).first << " ]" << std::endl;
      handle(o, (*i).second, pmf);
    }
  }
  template <typename O, typename C, typename F> void handleRMap_(const O* o, const C& c, F pmf) {
    for (typename C::const_iterator i = c.begin(); i != c.end(); ++i)  {
      const auto& cc = (*i).second;
      for (const auto& j : cc)   {
        (o->*pmf)(j);
      }
    }
  }

  std::string make_NCName(const std::string& in)   {
    std::string res = detail::str_replace(in, "/", "_");
    res = detail::str_replace(res, "#", "_");
    return res;
  }

  bool is_left_handed(const TGeoMatrix* m)   {
    const Double_t* r = m->GetRotationMatrix();
    if ( r )    {
      Double_t det =
        r[0]*r[4]*r[8] + r[3]*r[7]*r[2] + r[6]*r[1]*r[5] -
        r[2]*r[4]*r[6] - r[5]*r[7]*r[0] - r[8]*r[1]*r[3];
      return det < 0e0;
    }
    return false;
  }

  class G4UserRegionInformation : public G4VUserRegionInformation {
  public:
    Region region;
    double threshold;
    bool   storeSecondaries;
    G4UserRegionInformation()
      : threshold(0.0), storeSecondaries(false) {
    }
    virtual ~G4UserRegionInformation() {
    }
    virtual void Print() const {
      if (region.isValid())
        printout(DEBUG, "Region", "Name:%s", region.name());
    }
  };

  std::pair<double,double> g4PropertyConversion(int index)   {
#if G4VERSION_NUMBER >= 1040
    switch(index)  {
    case kRINDEX:                         return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kREFLECTIVITY:                   return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kREALRINDEX:                     return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kIMAGINARYRINDEX:                return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kEFFICIENCY:                     return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kTRANSMITTANCE:                  return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kSPECULARLOBECONSTANT:           return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kSPECULARSPIKECONSTANT:          return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kBACKSCATTERCONSTANT:            return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kGROUPVEL:                       return std::make_pair(CLHEP::keV/units::keV, (CLHEP::m/CLHEP::s)/(units::m/units::s));  // meter/second
    case kMIEHG:                          return std::make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
    case kRAYLEIGH:                       return std::make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);  // ??? says its a length
    case kWLSCOMPONENT:                   return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kWLSABSLENGTH:                   return std::make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
    case kABSLENGTH:                      return std::make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
#if G4VERSION_NUMBER >= 1100
    case kWLSCOMPONENT2:                  return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kWLSABSLENGTH2:                  return std::make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
    case kSCINTILLATIONCOMPONENT1:        return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
    case kSCINTILLATIONCOMPONENT2:        return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
    case kSCINTILLATIONCOMPONENT3:        return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
#else
    case kFASTCOMPONENT:                  return std::make_pair(CLHEP::keV/units::keV, 1.0);
    case kSLOWCOMPONENT:                  return std::make_pair(CLHEP::keV/units::keV, 1.0);
#endif
    case kPROTONSCINTILLATIONYIELD:       return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV); // Yields: 1/energy
    case kDEUTERONSCINTILLATIONYIELD:     return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
    case kTRITONSCINTILLATIONYIELD:       return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
    case kALPHASCINTILLATIONYIELD:        return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
    case kIONSCINTILLATIONYIELD:          return std::make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);