DDDividedPolycone.cc

Go to the documentation of this file.

#include "DetectorDescription/Parser/src/DDDividedPolycone.h"
#include "DetectorDescription/Core/interface/DDRotationMatrix.h"
#include "DetectorDescription/Core/interface/DDLogicalPart.h"
#include "DetectorDescription/Core/interface/DDMaterial.h"
#include "DetectorDescription/Core/interface/DDName.h"
#include "DetectorDescription/Core/interface/DDSolid.h"
#include "DetectorDescription/Core/interface/DDTransform.h"
#include "DataFormats/Math/interface/GeantUnits.h"
#include "DetectorDescription/Parser/src/DDDividedGeometryObject.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Utilities/interface/Exception.h"

#include <cstddef>
#include <iostream>
#include <string>
#include <utility>
#include <vector>

class DDCompactView;

using namespace geant_units::operators;

DDDividedPolyconeRho::DDDividedPolyconeRho(const DDDivision& div, DDCompactView* cpv)
    : DDDividedGeometryObject::DDDividedGeometryObject(div, cpv) {
  checkParametersValidity();
  DDPolycone msol = (DDPolycone)(div_.parent().solid());
  std::vector<double> localrMaxVec = msol.rMaxVec();
  std::vector<double> localrMinVec = msol.rMinVec();

  setType("DivisionPolyconeRho");

  // in DDD we only have ONE representation
  // in the case when rMinVec and rMaxVec
  // are empty rVec and zVec are the r and z std::vectors.

  if (divisionType_ == DivWIDTH) {
    compNDiv_ = calculateNDiv(localrMaxVec[0] - localrMinVec[0], div_.width(), div_.offset());
  } else if (divisionType_ == DivNDIV) {
    compWidth_ = calculateWidth(localrMaxVec[0] - localrMinVec[0], div_.nReplicas(), div_.offset());
  }
}

void DDDividedPolyconeRho::checkParametersValidity(void) {
  DDDividedGeometryObject::checkParametersValidity();

  DDPolycone msol = (DDPolycone)(div_.parent().solid());

  if (divisionType_ == DivNDIVandWIDTH || divisionType_ == DivWIDTH) {
    std::cout << "WARNING - "
              << "DDDividedPolyconeRho::checkParametersValidity()" << std::endl
              << "          Solid " << msol << std::endl
              << "          Division along r will be done with a width "
              << "different for each solid section." << std::endl
              << "          WIDTH will not be used !" << std::endl;
  }
  if (div_.offset() != 0.) {
    std::cout << "WARNING - "
              << "DDDividedPolyconeRho::checkParametersValidity()" << std::endl
              << "          Solid " << msol << std::endl
              << "          Division along  R will be done with a width "
              << "different for each solid section." << std::endl
              << "          OFFSET will not be used !" << std::endl;
  }
}

double DDDividedPolyconeRho::getMaxParameter(void) const {
  DDPolycone msol = (DDPolycone)(div_.parent().solid());
  std::vector<double> localrMaxVec = msol.rMaxVec();
  std::vector<double> localrMinVec = msol.rMinVec();

  return localrMaxVec[0] - localrMinVec[0];
}

DDRotation DDDividedPolyconeRho::makeDDRotation(const int copyNo) const {
  DDRotation myddrot;  // sets to identity.
  return myddrot;
}

DDTranslation DDDividedPolyconeRho::makeDDTranslation(const int copyNo) const {
  DDTranslation translation;
  return translation;
}

DDLogicalPart DDDividedPolyconeRho::makeDDLogicalPart(const int copyNo) const {
  DDSolid ddpolycone;
  DDMaterial usemat(div_.parent().material());

  DDPolycone msol = (DDPolycone)(div_.parent().solid());
  std::vector<double> localrMaxVec = msol.rMaxVec();
  std::vector<double> localrMinVec = msol.rMinVec();
  std::vector<double> localzVec = msol.zVec();

  int nZplanes = localzVec.size();

  std::vector<double> newrMinVec;
  std::vector<double> newrMaxVec;

  double width = 0.;
  for (int ii = 0; ii < nZplanes; ++ii) {
    width = calculateWidth(localrMaxVec[ii] - localrMinVec[ii], compNDiv_, div_.offset());
    newrMinVec.emplace_back(localrMinVec[ii] + div_.offset() + width * copyNo);
    newrMaxVec.emplace_back(localrMinVec[ii] + div_.offset() + width * (copyNo + 1));
  }
  DDName solname(div_.parent().ddname().name() + "_DIVCHILD" + std::to_string(copyNo), div_.parent().ddname().ns());

  ddpolycone = DDSolidFactory::polycone(solname, msol.startPhi(), msol.deltaPhi(), localzVec, newrMinVec, newrMaxVec);

  DDLogicalPart ddlp = DDLogicalPart(solname, usemat, ddpolycone);
  return ddlp;
}

DDDividedPolyconePhi::DDDividedPolyconePhi(const DDDivision& div, DDCompactView* cpv)
    : DDDividedGeometryObject::DDDividedGeometryObject(div, cpv) {
  checkParametersValidity();
  DDPolycone msol = (DDPolycone)(div_.parent().solid());
  setType("DivisionPolyconePhi");
  // this is the g4.  what do we keep??? I think it is deltaPhi
  // double deltaPhi = msol->GetEndPhi() - msol->GetStartPhi();
  if (divisionType_ == DivWIDTH) {
    //If you divide a tube of 360 degrees the offset displaces the starting angle, but you still fill the 360 degrees
    if (msol.deltaPhi() == 360._deg) {
      compNDiv_ = calculateNDiv(msol.deltaPhi(), div_.width(), 0.);
    } else {
      compNDiv_ = calculateNDiv(msol.deltaPhi(), div_.width(), div_.offset());
    }
  } else if (divisionType_ == DivNDIV) {
    if (msol.deltaPhi() == 360._deg) {
      compWidth_ = calculateWidth(msol.deltaPhi(), div_.nReplicas(), 0.);
    } else {
      compWidth_ = calculateWidth(msol.deltaPhi(), div_.nReplicas(), div_.offset());
    }
  }
}

void DDDividedPolyconePhi::checkParametersValidity(void) { DDDividedGeometryObject::checkParametersValidity(); }

double DDDividedPolyconePhi::getMaxParameter(void) const {
  DDPolycone msol = (DDPolycone)(div_.parent().solid());
  return msol.deltaPhi();
}

DDRotation DDDividedPolyconePhi::makeDDRotation(const int copyNo) const {
  DDRotation myddrot;  // sets to identity.
  double posi = (copyNo - 1) * compWidth_;
  // how to name the rotation??
  // i do not like this
  DDName ddrotname(div_.parent().ddname().name() + "_DIVCHILD_ROT" + std::to_string(copyNo),
                   div_.parent().ddname().ns());
  myddrot = DDrot(ddrotname, changeRotMatrix(posi));

  return myddrot;
}

DDTranslation DDDividedPolyconePhi::makeDDTranslation(const int copyNo) const {
  DDTranslation translation;
  return translation;
}

DDLogicalPart DDDividedPolyconePhi::makeDDLogicalPart(const int copyNo) const {
  DDSolid ddpolycone;
  DDMaterial usemat(div_.parent().material());

  DDPolycone msol = (DDPolycone)(div_.parent().solid());
  std::vector<double> localrMaxVec = msol.rMaxVec();
  std::vector<double> localrMinVec = msol.rMinVec();
  std::vector<double> localzVec = msol.zVec();

  DDName solname(div_.parent().ddname().name() + "_DIVCHILD", div_.parent().ddname().ns());
  DDSolid sol(solname);
  if (!sol.isDefined().second) {
    ddpolycone = DDSolidFactory::polycone(
        solname, msol.startPhi() + div_.offset(), compWidth_, localzVec, localrMinVec, localrMaxVec);
  }
  DDLogicalPart ddlp(solname);
  if (!ddlp.isDefined().second) {
    ddlp = DDLogicalPart(solname, usemat, ddpolycone);
  }

  return ddlp;
}

DDDividedPolyconeZ::DDDividedPolyconeZ(const DDDivision& div, DDCompactView* cpv)
    : DDDividedGeometryObject::DDDividedGeometryObject(div, cpv) {
  checkParametersValidity();
  DDPolycone msol = (DDPolycone)(div_.parent().solid());
  std::vector<double> localrMaxVec = msol.rMaxVec();
  std::vector<double> localrMinVec = msol.rMinVec();
  std::vector<double> localzVec = msol.zVec();

  if (divisionType_ == DivWIDTH) {
    compNDiv_ = calculateNDiv(localzVec[localzVec.size() - 1] - localzVec[0], div_.width(), div_.offset());
  } else if (divisionType_ == DivNDIV) {
    compWidth_ = calculateNDiv(localzVec[localzVec.size() - 1] - localzVec[0], div_.nReplicas(), div_.offset());
  }
}

void DDDividedPolyconeZ::checkParametersValidity(void) {
  DDDividedGeometryObject::checkParametersValidity();

  DDPolycone msol = (DDPolycone)(div_.parent().solid());
  std::vector<double> localzVec = msol.zVec();
  // CHANGE FROM G4 a polycone can be divided in Z by specifying
  // nReplicas IF they happen to coincide with the number of
  // z plans.
  size_t tempNDiv = div_.nReplicas();
  if (tempNDiv == 0)
    tempNDiv = calculateNDiv(localzVec[localzVec.size() - 1] - localzVec[0], div_.width(), div_.offset());
  if ((msol.zVec().size() - 1) != tempNDiv) {
    std::string s = "ERROR - DDDividedPolyconeZ::checkParametersValidity()";
    s += "\n\tDivision along Z will be done splitting in the defined";
    s += "\n\tz_planes, i.e, the number of division would be :";
    s += "\n\t" + std::to_string(msol.zVec().size() - 1);
    s += "\n\tinstead of " + std::to_string(tempNDiv) + " !\n";

    throw cms::Exception("DDException") << s;
  }
}

double DDDividedPolyconeZ::getMaxParameter(void) const {
  DDPolycone msol = (DDPolycone)(div_.parent().solid());
  std::vector<double> localzVec = msol.zVec();

  return (localzVec[localzVec.size() - 1] - localzVec[0]);
}

DDRotation DDDividedPolyconeZ::makeDDRotation(const int copyNo) const {
  DDRotation myddrot;  // sets to identity.
  return myddrot;
}

DDTranslation DDDividedPolyconeZ::makeDDTranslation(const int copyNo) const {
  DDTranslation translation;
  DDPolycone msol = (DDPolycone)(div_.parent().solid());
  std::vector<double> localzVec = msol.zVec();
  double posi = (localzVec[copyNo] + localzVec[copyNo + 1]) / 2;
  translation.SetZ(posi);
  return translation;
}

DDLogicalPart DDDividedPolyconeZ::makeDDLogicalPart(const int copyNo) const {
  DDSolid ddpolycone;
  DDMaterial usemat(div_.parent().material());

  DDPolycone msol = (DDPolycone)(div_.parent().solid());
  std::vector<double> localrMaxVec = msol.rMaxVec();
  std::vector<double> localrMinVec = msol.rMinVec();
  std::vector<double> localzVec = msol.zVec();

  DDName solname(div_.parent().ddname().name() + "_DIVCHILD" + std::to_string(copyNo), div_.parent().ddname().ns());
  ddpolycone = DDSolidFactory::cons(solname,
                                    compWidth_ / 2,
                                    localrMinVec[copyNo],
                                    localrMaxVec[copyNo],
                                    localrMinVec[copyNo + 1],
                                    localrMaxVec[copyNo + 1],
                                    msol.startPhi(),
                                    msol.deltaPhi());

  DDLogicalPart ddlp = DDLogicalPart(solname, usemat, ddpolycone);

  return ddlp;
}