TGeoFromDddService.cc

Go to the documentation of this file.

// -*- C++ -*-
//
// Package:     Geometry
// Class  :     TGeoFromDddService
//
// Implementation:
//     [Notes on implementation]
//
// Original Author:
//         Created:  Fri Jul  2 16:11:42 CEST 2010
//

// system include files

// user include files

#include "Fireworks/Geometry/interface/TGeoFromDddService.h"

#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/ActivityRegistry.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/ESTransientHandle.h"

#include "DetectorDescription/Core/interface/DDCompactView.h"
#include "DetectorDescription/Core/interface/DDSolid.h"
#include "DetectorDescription/Core/interface/DDMaterial.h"
#include "DataFormats/Math/interface/GraphWalker.h"

#include "Geometry/Records/interface/IdealGeometryRecord.h"

#include "TGeoManager.h"
#include "TGeoMatrix.h"
#include "TGeoCompositeShape.h"
#include "TGeoPcon.h"
#include "TGeoPgon.h"
#include "TGeoCone.h"
#include "TGeoBoolNode.h"
#include "TGeoTube.h"
#include "TGeoArb8.h"
#include "TGeoTrd2.h"
#include "TGeoXtru.h"

#include "Math/GenVector/RotationX.h"

#include "CLHEP/Units/GlobalSystemOfUnits.h"

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
TGeoFromDddService::TGeoFromDddService(const edm::ParameterSet& pset, edm::ActivityRegistry& ar)
    : m_level(pset.getUntrackedParameter<int>("level", 10)),
      m_verbose(pset.getUntrackedParameter<bool>("verbose", false)),
      m_eventSetup(nullptr),
      m_geoManager(nullptr) {
  ar.watchPostBeginRun(this, &TGeoFromDddService::postBeginRun);
  ar.watchPostEndRun(this, &TGeoFromDddService::postEndRun);
}

TGeoFromDddService::~TGeoFromDddService() {
  if (m_geoManager) {
    delete m_geoManager;
  }
}

//==============================================================================
// public member functions
//==============================================================================

void TGeoFromDddService::postBeginRun(const edm::Run&, const edm::EventSetup& es) {
  printf("TGeoFromDddService::postBeginRun\n");

  m_eventSetup = &es;
}

void TGeoFromDddService::postEndRun(const edm::Run&, const edm::EventSetup&) {
  printf("TGeoFromDddService::postEndRun\n");

  // Construction of geometry fails miserably on second attempt ...
  /*
   if (m_geoManager)
   {
      delete m_geoManager;
      m_geoManager = 0;
   }
   */
  m_eventSetup = nullptr;
}

TGeoManager* TGeoFromDddService::getGeoManager() {
  if (m_geoManager == nullptr) {
    if (m_eventSetup == nullptr)
      edm::LogError("TGeoFromDddService") << "getGeoManager() -- EventSetup not present.\n";
    else {
      m_geoManager = createManager(m_level);
      if (m_geoManager == nullptr)
        edm::LogError("TGeoFromDddService") << "getGeoManager() -- creation failed.\n";
    }
  }
  gGeoManager = m_geoManager;
  return m_geoManager;
}

//==============================================================================
// Local helpers
//==============================================================================

namespace {
  TGeoCombiTrans* createPlacement(const DDRotationMatrix& iRot, const DDTranslation& iTrans) {
    //  std::cout << "in createPlacement" << std::endl;
    double elements[9];
    iRot.GetComponents(elements);
    TGeoRotation r;
    r.SetMatrix(elements);

    TGeoTranslation t(iTrans.x() / cm, iTrans.y() / cm, iTrans.z() / cm);

    return new TGeoCombiTrans(t, r);
  }
}  // namespace

//==============================================================================
// private member functions
//==============================================================================

TGeoManager* TGeoFromDddService::createManager(int level) {
  using namespace edm;

  ESTransientHandle<DDCompactView> viewH;
  m_eventSetup->get<IdealGeometryRecord>().get(viewH);

  if (!viewH.isValid()) {
    return nullptr;
  }

  TGeoManager* geo_mgr = new TGeoManager("cmsGeo", "CMS Detector");
  // NOTE: the default constructor does not create the identity matrix
  if (gGeoIdentity == nullptr) {
    gGeoIdentity = new TGeoIdentity("Identity");
  }

  std::cout << "about to initialize the DDCompactView walker" << std::endl;
  auto walker = math::GraphWalker<DDLogicalPart, DDPosData*>(viewH->graph(), viewH->root());
  auto info = walker.current();

  // The top most item is actually the volume holding both the
  // geometry AND the magnetic field volumes!
  walker.firstChild();

  TGeoVolume* top = createVolume(info.first.name().fullname(), info.first.solid(), info.first.material());
  if (top == nullptr) {
    return nullptr;
  }

  geo_mgr->SetTopVolume(top);
  // ROOT chokes unless colors are assigned
  top->SetVisibility(kFALSE);
  top->SetLineColor(kBlue);

  std::vector<TGeoVolume*> parentStack;
  parentStack.push_back(top);

  if (!walker.firstChild()) {
    return nullptr;
  }

  do {
    auto info = walker.current();

    if (m_verbose) {
      for (unsigned int i = 0; i < parentStack.size(); ++i) {
        std::cout << " ";
      }
      std::cout << info.first.name() << " " << info.second->copyno() << " "
                << DDSolidShapesName::name(info.first.solid().shape()) << std::endl;
    }

    bool childAlreadyExists = (nullptr != nameToVolume_[info.first.name().fullname()]);
    TGeoVolume* child = createVolume(info.first.name().fullname(), info.first.solid(), info.first.material());
    if (nullptr != child && info.second != nullptr) {
      parentStack.back()->AddNode(
          child, info.second->copyno(), createPlacement(info.second->rotation(), info.second->translation()));
      child->SetLineColor(kBlue);
    } else {
      if (info.second == nullptr) {
        break;
      }
    }
    if (nullptr == child || childAlreadyExists || level == int(parentStack.size())) {
      if (nullptr != child) {
        child->SetLineColor(kRed);
      }
      //stop descending
      if (!walker.nextSibling()) {
        while (walker.parent()) {
          parentStack.pop_back();
          if (walker.nextSibling()) {
            break;
          }
        }
      }
    } else {
      if (walker.firstChild()) {
        parentStack.push_back(child);
      } else {
        if (!walker.nextSibling()) {
          while (walker.parent()) {
            parentStack.pop_back();
            if (walker.nextSibling()) {
              break;
            }
          }
        }
      }
    }
  } while (!parentStack.empty());

  geo_mgr->CloseGeometry();

  geo_mgr->DefaultColors();

  nameToShape_.clear();
  nameToVolume_.clear();
  nameToMaterial_.clear();
  nameToMedium_.clear();

  return geo_mgr;
}

TGeoShape* TGeoFromDddService::createShape(const std::string& iName, const DDSolid& iSolid) {
  TGeoShape* rSolid = nameToShape_[iName];
  if (rSolid == nullptr) {
    const std::vector<double>& params = iSolid.parameters();
    //      std::cout <<"  shape "<<iSolid<<std::endl;
    switch (iSolid.shape()) {
      case DDSolidShape::ddbox:
        rSolid = new TGeoBBox(iName.c_str(), params[0] / cm, params[1] / cm, params[2] / cm);
        break;
      case DDSolidShape::ddcons:
        rSolid = new TGeoConeSeg(iName.c_str(),
                                 params[0] / cm,
                                 params[1] / cm,
                                 params[2] / cm,
                                 params[3] / cm,
                                 params[4] / cm,
                                 params[5] / deg,
                                 params[6] / deg + params[5] / deg);
        break;
      case DDSolidShape::ddtubs:
        //Order in params is  zhalf,rIn,rOut,startPhi,deltaPhi
        rSolid = new TGeoTubeSeg(iName.c_str(),
                                 params[1] / cm,
                                 params[2] / cm,
                                 params[0] / cm,
                                 params[3] / deg,
                                 params[3] / deg + params[4] / deg);
        break;
      case DDSolidShape::ddcuttubs:
        //Order in params is  zhalf,rIn,rOut,startPhi,deltaPhi
        rSolid = new TGeoCtub(iName.c_str(),
                              params[1] / cm,
                              params[2] / cm,
                              params[0] / cm,
                              params[3] / deg,
                              params[3] / deg + params[4] / deg,
                              params[5],
                              params[6],
                              params[7],
                              params[8],
                              params[9],
                              params[10]);
        break;
      case DDSolidShape::ddtrap:
        rSolid = new TGeoTrap(iName.c_str(),
                              params[0] / cm,     //dz
                              params[1] / deg,    //theta
                              params[2] / deg,    //phi
                              params[3] / cm,     //dy1
                              params[4] / cm,     //dx1
                              params[5] / cm,     //dx2
                              params[6] / deg,    //alpha1
                              params[7] / cm,     //dy2
                              params[8] / cm,     //dx3
                              params[9] / cm,     //dx4
                              params[10] / deg);  //alpha2
        break;
      case DDSolidShape::ddpolycone_rrz:
        rSolid = new TGeoPcon(iName.c_str(), params[0] / deg, params[1] / deg, (params.size() - 2) / 3);
        {
          std::vector<double> temp(params.size() + 1);
          temp.reserve(params.size() + 1);
          temp[0] = params[0] / deg;
          temp[1] = params[1] / deg;
          temp[2] = (params.size() - 2) / 3;
          std::copy(params.begin() + 2, params.end(), temp.begin() + 3);
          for (std::vector<double>::iterator it = temp.begin() + 3; it != temp.end(); ++it) {
            *it /= cm;
          }
          rSolid->SetDimensions(&(*(temp.begin())));
        }
        break;
      case DDSolidShape::ddpolyhedra_rrz:
        rSolid = new TGeoPgon(
            iName.c_str(), params[1] / deg, params[2] / deg, static_cast<int>(params[0]), (params.size() - 3) / 3);
        {
          std::vector<double> temp(params.size() + 1);
          temp[0] = params[1] / deg;
          temp[1] = params[2] / deg;
          temp[2] = params[0];
          temp[3] = (params.size() - 3) / 3;
          std::copy(params.begin() + 3, params.end(), temp.begin() + 4);
          for (std::vector<double>::iterator it = temp.begin() + 4; it != temp.end(); ++it) {
            *it /= cm;
          }
          rSolid->SetDimensions(&(*(temp.begin())));
        }
        break;
      case DDSolidShape::ddextrudedpolygon: {
        DDExtrudedPolygon extrPgon(iSolid);
        std::vector<double> x = extrPgon.xVec();
        std::transform(x.begin(), x.end(), x.begin(), [](double d) { return d / cm; });
        std::vector<double> y = extrPgon.yVec();
        std::transform(y.begin(), y.end(), y.begin(), [](double d) { return d / cm; });
        std::vector<double> z = extrPgon.zVec();
        std::vector<double> zx = extrPgon.zxVec();
        std::vector<double> zy = extrPgon.zyVec();
        std::vector<double> zscale = extrPgon.zscaleVec();

        TGeoXtru* mySolid = new TGeoXtru(z.size());
        mySolid->DefinePolygon(x.size(), &(*x.begin()), &(*y.begin()));
        for (size_t i = 0; i < params[0]; ++i) {
          mySolid->DefineSection(i, z[i] / cm, zx[i] / cm, zy[i] / cm, zscale[i]);
        }

        rSolid = mySolid;
      } break;
      case DDSolidShape::ddpseudotrap: {
        //implementation taken from SimG4Core/Geometry/src/DDG4SolidConverter.cc
        const static DDRotationMatrix s_rot(ROOT::Math::RotationX(90. * deg));
        DDPseudoTrap pt(iSolid);
        assert(pt.radius() < 0);
        double x = 0;
        double r = fabs(pt.radius());
        if (pt.atMinusZ()) {
          x = pt.x1();
        } else {
          x = pt.x2();
        }
        double openingAngle = 2.0 * asin(x / r);
        double h = pt.y1() < pt.y2() ? pt.y2() : pt.y1();
        h += h / 20.;
        double displacement = 0;
        double startPhi = 0;
        double delta = sqrt((r + x) * (r - x));
        if (pt.atMinusZ()) {
          displacement = -pt.halfZ() - delta;
          startPhi = 270. - openingAngle / deg / 2.0;
        } else {
          displacement = pt.halfZ() + delta;
          startPhi = 90. - openingAngle / deg / 2.;
        }
        std::unique_ptr<TGeoShape> trap(new TGeoTrd2(
            pt.name().name().c_str(), pt.x1() / cm, pt.x2() / cm, pt.y1() / cm, pt.y2() / cm, pt.halfZ() / cm));
        std::unique_ptr<TGeoShape> tubs(
            new TGeoTubeSeg(pt.name().name().c_str(), 0., r / cm, h / cm, startPhi, openingAngle));
        TGeoSubtraction* sub = new TGeoSubtraction(
            trap.release(), tubs.release(), createPlacement(s_rot, DDTranslation(0., 0., displacement)));
        rSolid = new TGeoCompositeShape(iName.c_str(), sub);

        break;
      }
      case DDSolidShape::ddsubtraction: {
        DDBooleanSolid boolSolid(iSolid);
        if (!boolSolid) {
          throw cms::Exception("GeomConvert") << "conversion to DDBooleanSolid failed";
        }

        std::unique_ptr<TGeoShape> left(createShape(boolSolid.solidA().name().fullname(), boolSolid.solidA()));
        std::unique_ptr<TGeoShape> right(createShape(boolSolid.solidB().name().fullname(), boolSolid.solidB()));
        if (nullptr != left.get() && nullptr != right.get()) {
          TGeoSubtraction* sub =
              new TGeoSubtraction(left.release(),
                                  right.release(),
                                  gGeoIdentity,
                                  createPlacement(boolSolid.rotation().matrix(), boolSolid.translation()));
          rSolid = new TGeoCompositeShape(iName.c_str(), sub);
        }
        break;
      }
      case DDSolidShape::ddunion: {
        DDBooleanSolid boolSolid(iSolid);
        if (!boolSolid) {
          throw cms::Exception("GeomConvert") << "conversion to DDBooleanSolid failed";
        }

        std::unique_ptr<TGeoShape> left(createShape(boolSolid.solidA().name().fullname(), boolSolid.solidA()));
        std::unique_ptr<TGeoShape> right(createShape(boolSolid.solidB().name().fullname(), boolSolid.solidB()));
        //DEBUGGING
        //break;
        if (nullptr != left.get() && nullptr != right.get()) {
          TGeoUnion* boolS = new TGeoUnion(left.release(),
                                           right.release(),
                                           gGeoIdentity,
                                           createPlacement(boolSolid.rotation().matrix(), boolSolid.translation()));
          rSolid = new TGeoCompositeShape(iName.c_str(), boolS);
        }
        break;
      }
      case DDSolidShape::ddintersection: {
        DDBooleanSolid boolSolid(iSolid);
        if (!boolSolid) {
          throw cms::Exception("GeomConvert") << "conversion to DDBooleanSolid failed";
        }

        std::unique_ptr<TGeoShape> left(createShape(boolSolid.solidA().name().fullname(), boolSolid.solidA()));
        std::unique_ptr<TGeoShape> right(createShape(boolSolid.solidB().name().fullname(), boolSolid.solidB()));
        if (nullptr != left.get() && nullptr != right.get()) {
          TGeoIntersection* boolS =
              new TGeoIntersection(left.release(),
                                   right.release(),
                                   gGeoIdentity,
                                   createPlacement(boolSolid.rotation().matrix(), boolSolid.translation()));
          rSolid = new TGeoCompositeShape(iName.c_str(), boolS);
        }
        break;
      }
      default:
        break;
    }
    nameToShape_[iName] = rSolid;
  }
  if (rSolid == nullptr) {
    std::cerr << "COULD NOT MAKE " << iName << std::endl;
  }
  return rSolid;
}

TGeoVolume* TGeoFromDddService::createVolume(const std::string& iName,
                                             const DDSolid& iSolid,
                                             const DDMaterial& iMaterial) {
  TGeoVolume* v = nameToVolume_[iName];
  if (v == nullptr) {
    TGeoShape* solid = createShape(iSolid.name().fullname(), iSolid);
    std::string mat_name = iMaterial.name().fullname();
    TGeoMedium* geo_med = nameToMedium_[mat_name];
    if (geo_med == nullptr) {
      TGeoMaterial* geo_mat = createMaterial(iMaterial);
      geo_med = new TGeoMedium(mat_name.c_str(), 0, geo_mat);
      nameToMedium_[mat_name] = geo_med;
    }
    if (solid) {
      v = new TGeoVolume(iName.c_str(), solid, geo_med);
    }
    nameToVolume_[iName] = v;
  }
  return v;
}

TGeoMaterial* TGeoFromDddService::createMaterial(const DDMaterial& iMaterial) {
  std::string mat_name = iMaterial.name().fullname();
  TGeoMaterial* mat = nameToMaterial_[mat_name];

  if (mat == nullptr) {
    if (iMaterial.noOfConstituents() > 0) {
      TGeoMixture* mix = new TGeoMixture(mat_name.c_str(), iMaterial.noOfConstituents(), iMaterial.density() * cm3 / g);
      for (int i = 0; i < iMaterial.noOfConstituents(); ++i) {
        mix->AddElement(createMaterial(iMaterial.constituent(i).first), iMaterial.constituent(i).second);
      }
      mat = mix;
    } else {
      mat = new TGeoMaterial(mat_name.c_str(), iMaterial.a() * mole / g, iMaterial.z(), iMaterial.density() * cm3 / g);
    }
    nameToMaterial_[mat_name] = mat;
  }

  return mat;
}

//
// const member functions
//

//
// static member functions
//