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/Detector.h"
#include "DD4hep/Plugins.h"
#include "DD4hep/Shapes.h"
#include "DD4hep/Volumes.h"
#include "DD4hep/Printout.h"
#include "DD4hep/DD4hepUnits.h"
#include "DD4hep/PropertyTable.h"
#include "DD4hep/detail/ShapesInterna.h"
#include "DD4hep/detail/ObjectsInterna.h"
#include "DD4hep/detail/DetectorInterna.h"

#include "DDG4/Geant4Field.h"
#include "DDG4/Geant4Converter.h"
#include "DDG4/Geant4UserLimits.h"
#include "Geant4ShapeConverter.h"

// ROOT includes
#include "TMath.h"
#include "TROOT.h"
//#include "TColor.h"
//#include "TGeoManager.h"
#include "TGeoBoolNode.h"

// Geant4 include files
#include "G4Version.hh"
#include "G4VisAttributes.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 "G4Element.hh"
#include "G4Isotope.hh"
#include "G4Material.hh"
#include "G4UserLimits.hh"
#include "G4FieldManager.hh"
#include "G4LogicalVolume.hh"
#include "G4ReflectionFactory.hh"
#include "G4OpticalSurface.hh"
#include "G4LogicalSkinSurface.hh"
#include "G4ElectroMagneticField.hh"
#include "G4LogicalBorderSurface.hh"
#include "G4MaterialPropertiesTable.hh"
#if G4VERSION_NUMBER >= 1040
#include "G4MaterialPropertiesIndex.hh"
#endif
#if G4VERSION_NUMBER >= 1030
#include "G4ScaledSolid.hh"
#endif
#include "CLHEP/Units/SystemOfUnits.h"

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

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

#include "DDG4/Geant4AssemblyVolume.h"
#include "DD4hep/DetectorTools.h"

static const double CM_2_MM = (CLHEP::centimeter/dd4hep::centimeter);

namespace {
  static string indent = "";
  static Double_t s_identity_rot[] = { 1., 0., 0., 0., 1., 0., 0., 0., 1. };
  struct MyTransform3D : public G4Transform3D {
    MyTransform3D(double XX, double XY, double XZ, double DX, double YX, double YY, double YZ, double DY, double ZX, double ZY,
                  double ZZ, double DZ)
      : G4Transform3D(XX, XY, XZ, DX, YX, YY, YZ, DY, ZX, ZY, ZZ, DZ) {
    }
    MyTransform3D(const double* t, const double* r = s_identity_rot)
      : G4Transform3D(r[0],r[1],r[2],t[0]*CM_2_MM,r[3],r[4],r[5],t[1]*CM_2_MM,r[6],r[7],r[8],t[2]*CM_2_MM)  {
    }
    MyTransform3D(Transform3D&& copy) : Transform3D(copy) {}
  };

  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());
    }
  };


  pair<double,double> g4PropertyConversion(int index)   {
#if G4VERSION_NUMBER >= 1040
    switch(index)  {
    case kRINDEX:                         return make_pair(CLHEP::keV/units::keV, 1.0);
    case kREFLECTIVITY:                   return make_pair(CLHEP::keV/units::keV, 1.0);
    case kREALRINDEX:                     return make_pair(CLHEP::keV/units::keV, 1.0);
    case kIMAGINARYRINDEX:                return make_pair(CLHEP::keV/units::keV, 1.0);
    case kEFFICIENCY:                     return make_pair(CLHEP::keV/units::keV, 1.0);
    case kTRANSMITTANCE:                  return make_pair(CLHEP::keV/units::keV, 1.0);
    case kSPECULARLOBECONSTANT:           return make_pair(CLHEP::keV/units::keV, 1.0);
    case kSPECULARSPIKECONSTANT:          return make_pair(CLHEP::keV/units::keV, 1.0);
    case kBACKSCATTERCONSTANT:            return make_pair(CLHEP::keV/units::keV, 1.0);
    case kGROUPVEL:                       return make_pair(CLHEP::keV/units::keV, (CLHEP::m/CLHEP::s)/(units::m/units::s));  // meter/second
    case kMIEHG:                          return make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
    case kRAYLEIGH:                       return make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);  // ??? says its a length
    case kWLSCOMPONENT:                   return make_pair(CLHEP::keV/units::keV, 1.0);
    case kWLSABSLENGTH:                   return make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
    case kABSLENGTH:                      return make_pair(CLHEP::keV/units::keV, CLHEP::m/units::m);
    case kFASTCOMPONENT:                  return make_pair(CLHEP::keV/units::keV, 1.0);
    case kSLOWCOMPONENT:                  return make_pair(CLHEP::keV/units::keV, 1.0);
    case kPROTONSCINTILLATIONYIELD:       return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV); // Yields: 1/energy
    case kDEUTERONSCINTILLATIONYIELD:     return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
    case kTRITONSCINTILLATIONYIELD:       return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
    case kALPHASCINTILLATIONYIELD:        return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
    case kIONSCINTILLATIONYIELD:          return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
    case kELECTRONSCINTILLATIONYIELD:     return make_pair(CLHEP::keV/units::keV, units::keV/CLHEP::keV);
    default:
      break;
    }
    printout(FATAL,"Geant4Converter", "+++ Cannot convert material property with index: %d", index);
#else
    printout(FATAL,"Geant4Converter", "+++ Cannot convert material property with index: %d [Need Geant4 > 10.03]", index);
#endif
    return make_pair(0e0,0e0);
  }

  double g4ConstPropertyConversion(int index)   {
#if G4VERSION_NUMBER >= 1040
    switch(index)   {
    case kSURFACEROUGHNESS:            return CLHEP::m/units::m;                             // Length
    case kISOTHERMAL_COMPRESSIBILITY:  return (CLHEP::m3/CLHEP::keV)/(units::m3/CLHEP::keV); // Volume/Energy
    case kRS_SCALE_FACTOR:             return 1.0;  // ??
    case kWLSMEANNUMBERPHOTONS:        return 1.0;  // ??
    case kWLSTIMECONSTANT:             return CLHEP::second/units::second;                   // Time
    case kMIEHG_FORWARD:               return 1.0;
    case kMIEHG_BACKWARD:              return 1.0;
    case kMIEHG_FORWARD_RATIO:         return 1.0;
    case kSCINTILLATIONYIELD:          return units::keV/CLHEP::keV;                         // Energy
    case kRESOLUTIONSCALE:             return 1.0;
    case kFASTTIMECONSTANT:            return CLHEP::second/units::second;                   // Time
    case kFASTSCINTILLATIONRISETIME:   return CLHEP::second/units::second;                   // Time
    case kSLOWTIMECONSTANT:            return CLHEP::second/units::second;                   // Time
    case kSLOWSCINTILLATIONRISETIME:   return CLHEP::second/units::second;                   // Time
    case kYIELDRATIO:                  return 1.0;
    case kFERMIPOT:                    return CLHEP::keV/units::keV;                         // Energy
    case kDIFFUSION:                   return 1.0;
    case kSPINFLIP:                    return 1.0;
    case kLOSS:                        return 1.0;  // ??
    case kLOSSCS:                      return CLHEP::barn/units::barn;  // ??
    case kABSCS:                       return CLHEP::barn/units::barn;  // ??
    case kSCATCS:                      return CLHEP::barn/units::barn;  // ??
    case kMR_NBTHETA:                  return 1.0;
    case kMR_NBE:                      return 1.0;
    case kMR_RRMS:                     return 1.0;  // ??
    case kMR_CORRLEN:                  return CLHEP::m/units::m;                             // Length
    case kMR_THETAMIN:                 return 1.0;
    case kMR_THETAMAX:                 return 1.0;
    case kMR_EMIN:                     return CLHEP::keV/units::keV;                         // Energy