Parser/src/DDDividedPolyhedra.cc

0001 #include "DetectorDescription/Parser/src/DDDividedPolyhedra.h"
0002 #include "DetectorDescription/Core/interface/DDRotationMatrix.h"
0003 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
0004 #include "DetectorDescription/Core/interface/DDMaterial.h"
0005 #include "DetectorDescription/Core/interface/DDName.h"
0006 #include "DetectorDescription/Core/interface/DDSolid.h"
0007 #include "DetectorDescription/Core/interface/DDTransform.h"
0008 #include "DataFormats/Math/interface/GeantUnits.h"
0009 #include "DetectorDescription/Parser/src/DDDividedGeometryObject.h"
0010 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0011 #include "FWCore/Utilities/interface/Exception.h"
0012
0013 #include <cstddef>
0014 #include <iostream>
0015 #include <string>
0016 #include <utility>
0017 #include <vector>
0018
0019 class DDCompactView;
0020
0021 using namespace geant_units::operators;
0022
0023 DDDividedPolyhedraRho::DDDividedPolyhedraRho(const DDDivision& div, DDCompactView* cpv)
0024     : DDDividedGeometryObject(div, cpv) {
0025   checkParametersValidity();
0026   setType("DivisionPolyhedraRho");
0027
0028   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0029
0030   if (divisionType_ == DivWIDTH) {
0031     compNDiv_ = calculateNDiv(msol.rMaxVec()[0] - msol.rMinVec()[0], div_.width(), div_.offset());
0032   } else if (divisionType_ == DivNDIV) {
0033     compWidth_ = calculateWidth(msol.rMaxVec()[0] - msol.rMinVec()[0], div_.nReplicas(), div_.offset());
0034   }
0035 }
0036
0037 void DDDividedPolyhedraRho::checkParametersValidity(void) {
0038   DDDividedGeometryObject::checkParametersValidity();
0039
0040   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0041
0042   if (divisionType_ == DivNDIVandWIDTH || divisionType_ == DivWIDTH) {
0043     std::cout << "WARNING - "
0044               << "DDDividedPolyhedraRho::checkParametersValidity()" << std::endl
0045               << "          Solid " << msol << std::endl
0046               << "          Division along R will be done with a width "
0047               << "different for each solid section." << std::endl
0048               << "          WIDTH will not be used !" << std::endl;
0049   }
0050   if (div_.offset() != 0.) {
0051     std::cout << "WARNING - "
0052               << "DDDividedPolyhedraRho::checkParametersValidity()" << std::endl
0053               << "          Solid " << msol << std::endl
0054               << "          Division along  R will be done with a width "
0055               << "different for each solid section." << std::endl
0056               << "          OFFSET will not be used !" << std::endl;
0057   }
0058 }
0059
0060 double DDDividedPolyhedraRho::getMaxParameter(void) const {
0061   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0062   return msol.rMaxVec()[0] - msol.rMinVec()[0];
0063 }
0064
0065 DDTranslation DDDividedPolyhedraRho::makeDDTranslation(const int copyNo) const { return DDTranslation(); }
0066
0067 DDRotation DDDividedPolyhedraRho::makeDDRotation(const int copyNo) const { return DDRotation(); }
0068
0069 DDLogicalPart DDDividedPolyhedraRho::makeDDLogicalPart(const int copyNo) const {
0070   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0071   DDMaterial usemat = div_.parent().material();
0072
0073   std::vector<double> localrMaxVec = msol.rMaxVec();
0074   std::vector<double> localrMinVec = msol.rMinVec();
0075   std::vector<double> localzVec = msol.zVec();
0076   std::vector<double> newrMinVec;
0077   std::vector<double> newrMaxVec;
0078   int nZplanes = localzVec.size();
0079
0080   double width = 0.;
0081   for (int ii = 0; ii < nZplanes; ++ii) {
0082     //     width = CalculateWidth( origparamMother->Rmax[ii]
0083     //                           - origparamMother->Rmin[ii], compNDiv_, foffset );
0084     //     origparam.Rmin[ii] = origparamMother->Rmin[ii]+foffset+width*copyNo;
0085     //     origparam.Rmax[ii] = origparamMother->Rmin[ii]+foffset+width*(copyNo+1);
0086     width = calculateWidth(localrMaxVec[ii] - localrMinVec[ii], compNDiv_, div_.offset());
0087     newrMinVec[ii] = localrMinVec[ii] + div_.offset() + width * copyNo;
0088     newrMaxVec[ii] = localrMaxVec[ii] + div_.offset() + width * (copyNo + 1);
0089   }
0090
0091   //   phedra.SetOriginalParameters(&origparam); // copy values & transfer pointers
0092   //   phedra.Reset();                           // reset to new solid parameters
0093
0094   DDName solname(div_.parent().ddname().name() + "_DIVCHILD" + std::to_string(copyNo), div_.parent().ddname().ns());
0095
0096   DDSolid dsol = DDSolidFactory::polyhedra(
0097       solname, msol.sides(), msol.startPhi(), msol.deltaPhi(), localzVec, newrMinVec, newrMaxVec);
0098   DDLogicalPart ddlp = DDLogicalPart(solname, usemat, dsol);
0099   return ddlp;
0100 }
0101
0102 DDDividedPolyhedraPhi::DDDividedPolyhedraPhi(const DDDivision& div, DDCompactView* cpv)
0103     : DDDividedGeometryObject(div, cpv) {
0104   checkParametersValidity();
0105   setType("DivisionPolyhedraPhi");
0106
0107   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0108
0109   if (divisionType_ == DivWIDTH) {
0110     //If you divide a tube of 360 degrees the offset displaces the starting angle, but you still fill the 360 degrees
0111     if (msol.deltaPhi() == 360._deg) {
0112       compNDiv_ = calculateNDiv(msol.deltaPhi(), div_.width(), 0.);
0113     } else {
0114       compNDiv_ = calculateNDiv(msol.deltaPhi(), div_.width(), div_.offset());
0115     }
0116   } else if (divisionType_ == DivNDIV) {
0117     if (msol.deltaPhi() == 360._deg) {
0118       compWidth_ = calculateWidth(msol.deltaPhi(), div_.nReplicas(), 0.);
0119     } else {
0120       // original line looks wrong!
0121       compWidth_ = calculateWidth(msol.deltaPhi(), div_.nReplicas(), div_.offset());
0122     }
0123   }
0124 }
0125
0126 double DDDividedPolyhedraPhi::getMaxParameter(void) const {
0127   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0128   return msol.deltaPhi();  //msol->GetEndPhi() - msol->GetStartPhi();
0129 }
0130
0131 void DDDividedPolyhedraPhi::checkParametersValidity(void) {
0132   DDDividedGeometryObject::checkParametersValidity();
0133
0134   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0135
0136   if (divisionType_ == DivNDIVandWIDTH || divisionType_ == DivWIDTH) {
0137     std::cout << "WARNING - "
0138               << "DDDividedPolyhedraPhi::checkParametersValidity()" << std::endl
0139               << "          Solid " << msol << std::endl
0140               << "          Division along PHI will be done splitting "
0141               << "in the defined numSide." << std::endl
0142               << "          WIDTH will not be used !" << std::endl;
0143   }
0144   if (div_.offset() != 0.) {
0145     std::cout << "WARNING - "
0146               << "DDDividedPolyhedraPhi::checkParametersValidity()" << std::endl
0147               << "          Solid " << msol << std::endl
0148               << "          Division along PHI will be done splitting "
0149               << "in the defined numSide." << std::endl
0150               << "          OFFSET will not be used !" << std::endl;
0151   }
0152
0153   if (msol.sides() != compNDiv_) {
0154     std::cout << "ERROR - "
0155               << "DDDividedPolyhedraPhi::checkParametersValidity()" << std::endl
0156               << "        Division along PHI will be done splitting in the defined" << std::endl
0157               << "        numSide, i.e, the number of division would be :"
0158               << "        " << msol.sides() << " instead of " << compNDiv_ << " !" << std::endl;
0159     std::string s = "DDDividedPolyhedraPhi::checkParametersValidity() Not supported configuration.";
0160     throw cms::Exception("DDException") << s;
0161   }
0162 }
0163
0164 DDTranslation DDDividedPolyhedraPhi::makeDDTranslation(const int copyNo) const { return DDTranslation(); }
0165
0166 DDRotation DDDividedPolyhedraPhi::makeDDRotation(const int copyNo) const {
0167   DDRotation myddrot;  // sets to identity.
0168   double posi = (copyNo - 1) * compWidth_;
0169   DDName ddrotname(div_.parent().ddname().name() + "_DIVCHILD_ROT" + std::to_string(copyNo),
0170                    div_.parent().ddname().ns());
0171   myddrot = DDrot(ddrotname, changeRotMatrix(posi));
0172
0173   return myddrot;
0174 }
0175
0176 DDLogicalPart DDDividedPolyhedraPhi::makeDDLogicalPart(const int copyNo) const {
0177   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0178   DDMaterial usemat = div_.parent().material();
0179
0180   DDName solname(div_.parent().ddname().name() + "_DIVCHILD", div_.parent().ddname().ns());
0181   DDSolid dsol(solname);
0182   if (!dsol.isDefined().second) {
0183     dsol = DDSolidFactory::polyhedra(
0184         solname, msol.sides(), msol.startPhi() + div_.offset(), compWidth_, msol.zVec(), msol.rMinVec(), msol.rMaxVec());
0185   }
0186   DDLogicalPart ddlp(solname);
0187   if (!ddlp.isDefined().second)
0188     DDLogicalPart ddlp2 = DDLogicalPart(solname, usemat, dsol);
0189   return ddlp;
0190 }
0191
0192 DDDividedPolyhedraZ::DDDividedPolyhedraZ(const DDDivision& div, DDCompactView* cpv)
0193     : DDDividedGeometryObject(div, cpv) {
0194   checkParametersValidity();
0195   setType("DivisionPolyhedraZ");
0196
0197   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0198
0199   std::vector<double> zvec = msol.zVec();
0200
0201   if (divisionType_ == DivWIDTH) {
0202     compNDiv_ = calculateNDiv(zvec[zvec.size() - 1] - zvec[0], div_.width(), div_.offset());
0203   } else if (divisionType_ == DivNDIV) {
0204     compWidth_ = calculateWidth(zvec[zvec.size() - 1] - zvec[0], div_.nReplicas(), div_.offset());
0205   }
0206 }
0207
0208 double DDDividedPolyhedraZ::getMaxParameter(void) const {
0209   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0210
0211   std::vector<double> zvec = msol.zVec();
0212   return (zvec[zvec.size() - 1] - zvec[0]);
0213 }
0214
0215 void DDDividedPolyhedraZ::checkParametersValidity(void) {
0216   DDDividedGeometryObject::checkParametersValidity();
0217
0218   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0219
0220   if (divisionType_ == DivNDIVandWIDTH || divisionType_ == DivWIDTH) {
0221     std::cout << "WARNING - "
0222               << "DDDividedPolyhedraZ::checkParametersValidity()" << std::endl
0223               << "          Solid " << msol << std::endl
0224               << "          Division along Z will be done splitting "
0225               << "in the defined z_planes." << std::endl
0226               << "          WIDTH will not be used !" << std::endl;
0227   }
0228
0229   if (div_.offset() != 0.) {
0230     std::cout << "WARNING - "
0231               << "DDDividedPolyhedraZ::checkParametersValidity()" << std::endl
0232               << "          Solid " << msol << std::endl
0233               << "          Division along Z will be done splitting "
0234               << "in the defined z_planes." << std::endl
0235               << "          OFFSET will not be used !" << std::endl;
0236   }
0237
0238   std::vector<double> zvec = msol.zVec();
0239
0240   if (zvec.size() - 1 != size_t(compNDiv_)) {
0241     std::cout << "ERROR - "
0242               << "DDDividedPolyhedraZ::checkParametersValidity()" << std::endl
0243               << "        Division along Z can only be done by splitting in the defined" << std::endl
0244               << "        z_planes, i.e, the number of division would be :"
0245               << "        " << zvec.size() - 1 << " instead of " << compNDiv_ << " !" << std::endl;
0246     std::string s = "DDDividedPolyhedraZ::checkParametersValidity()";
0247     s += "Illegal Construct. Not a supported configuration.";
0248     throw cms::Exception("DDException") << s;
0249   }
0250 }
0251
0252 DDTranslation DDDividedPolyhedraZ::makeDDTranslation(const int copyNo) const {
0253   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0254   std::vector<double> zvec = msol.zVec();
0255
0256   //----- set translation: along Z axis
0257   double posi = (zvec[copyNo] + zvec[copyNo + 1]) / 2;
0258
0259   DDTranslation tr(0, 0, posi);
0260   //----- calculate rotation matrix: unit
0261   return tr;
0262 }
0263
0264 DDRotation DDDividedPolyhedraZ::makeDDRotation(const int copyNo) const { return DDRotation(); }
0265
0266 DDLogicalPart DDDividedPolyhedraZ::makeDDLogicalPart(const int copyNo) const {
0267   // only for mother number of planes = 2!!
0268   // mec: what?  why?  comment above and = 2 below straight from G4 impl.
0269   DDPolyhedra msol = (DDPolyhedra)(div_.parent().solid());
0270   DDMaterial usemat = div_.parent().material();
0271
0272   std::vector<double> zvec = msol.zVec();
0273   std::vector<double> rminvec = msol.rMinVec();
0274   std::vector<double> rmaxvec = msol.rMaxVec();
0275
0276   double posi = (zvec[copyNo] + zvec[copyNo + 1]) / 2.0;
0277
0278   DDName solname(div_.parent().ddname().name() + "_DIVCHILD" + std::to_string(copyNo), div_.parent().ddname().ns());
0279   std::vector<double> newRmin, newRmax, newZ;
0280   newZ.emplace_back(zvec[copyNo] - posi);
0281   newZ.emplace_back(zvec[copyNo + 1] - posi);
0282   newRmin.emplace_back(rminvec[copyNo]);
0283   newRmin.emplace_back(rminvec[copyNo + 1]);
0284   newRmax.emplace_back(rmaxvec[copyNo]);
0285   newRmax.emplace_back(rmaxvec[copyNo + 1]);
0286
0287   DDSolid dsol =
0288       DDSolidFactory::polyhedra(solname, msol.sides(), msol.startPhi(), msol.deltaPhi(), newZ, newRmin, newRmax);
0289   DDLogicalPart lp(solname, usemat, dsol);
0290   return lp;
0291 }