Parser/src/DDDividedPolycone.cc

0001 #include "DetectorDescription/Parser/src/DDDividedPolycone.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 DDDividedPolyconeRho::DDDividedPolyconeRho(const DDDivision& div, DDCompactView* cpv)
0024     : DDDividedGeometryObject::DDDividedGeometryObject(div, cpv) {
0025   checkParametersValidity();
0026   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0027   std::vector<double> localrMaxVec = msol.rMaxVec();
0028   std::vector<double> localrMinVec = msol.rMinVec();
0029
0030   setType("DivisionPolyconeRho");
0031
0032   // in DDD we only have ONE representation
0033   // in the case when rMinVec and rMaxVec
0034   // are empty rVec and zVec are the r and z std::vectors.
0035
0036   if (divisionType_ == DivWIDTH) {
0037     compNDiv_ = calculateNDiv(localrMaxVec[0] - localrMinVec[0], div_.width(), div_.offset());
0038   } else if (divisionType_ == DivNDIV) {
0039     compWidth_ = calculateWidth(localrMaxVec[0] - localrMinVec[0], div_.nReplicas(), div_.offset());
0040   }
0041 }
0042
0043 void DDDividedPolyconeRho::checkParametersValidity(void) {
0044   DDDividedGeometryObject::checkParametersValidity();
0045
0046   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0047
0048   if (divisionType_ == DivNDIVandWIDTH || divisionType_ == DivWIDTH) {
0049     std::cout << "WARNING - "
0050               << "DDDividedPolyconeRho::checkParametersValidity()" << std::endl
0051               << "          Solid " << msol << std::endl
0052               << "          Division along r will be done with a width "
0053               << "different for each solid section." << std::endl
0054               << "          WIDTH will not be used !" << std::endl;
0055   }
0056   if (div_.offset() != 0.) {
0057     std::cout << "WARNING - "
0058               << "DDDividedPolyconeRho::checkParametersValidity()" << std::endl
0059               << "          Solid " << msol << std::endl
0060               << "          Division along  R will be done with a width "
0061               << "different for each solid section." << std::endl
0062               << "          OFFSET will not be used !" << std::endl;
0063   }
0064 }
0065
0066 double DDDividedPolyconeRho::getMaxParameter(void) const {
0067   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0068   std::vector<double> localrMaxVec = msol.rMaxVec();
0069   std::vector<double> localrMinVec = msol.rMinVec();
0070
0071   return localrMaxVec[0] - localrMinVec[0];
0072 }
0073
0074 DDRotation DDDividedPolyconeRho::makeDDRotation(const int copyNo) const {
0075   DDRotation myddrot;  // sets to identity.
0076   return myddrot;
0077 }
0078
0079 DDTranslation DDDividedPolyconeRho::makeDDTranslation(const int copyNo) const {
0080   DDTranslation translation;
0081   return translation;
0082 }
0083
0084 DDLogicalPart DDDividedPolyconeRho::makeDDLogicalPart(const int copyNo) const {
0085   DDSolid ddpolycone;
0086   DDMaterial usemat(div_.parent().material());
0087
0088   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0089   std::vector<double> localrMaxVec = msol.rMaxVec();
0090   std::vector<double> localrMinVec = msol.rMinVec();
0091   std::vector<double> localzVec = msol.zVec();
0092
0093   int nZplanes = localzVec.size();
0094
0095   std::vector<double> newrMinVec;
0096   std::vector<double> newrMaxVec;
0097
0098   double width = 0.;
0099   for (int ii = 0; ii < nZplanes; ++ii) {
0100     width = calculateWidth(localrMaxVec[ii] - localrMinVec[ii], compNDiv_, div_.offset());
0101     newrMinVec.emplace_back(localrMinVec[ii] + div_.offset() + width * copyNo);
0102     newrMaxVec.emplace_back(localrMinVec[ii] + div_.offset() + width * (copyNo + 1));
0103   }
0104   DDName solname(div_.parent().ddname().name() + "_DIVCHILD" + std::to_string(copyNo), div_.parent().ddname().ns());
0105
0106   ddpolycone = DDSolidFactory::polycone(solname, msol.startPhi(), msol.deltaPhi(), localzVec, newrMinVec, newrMaxVec);
0107
0108   DDLogicalPart ddlp = DDLogicalPart(solname, usemat, ddpolycone);
0109   return ddlp;
0110 }
0111
0112 DDDividedPolyconePhi::DDDividedPolyconePhi(const DDDivision& div, DDCompactView* cpv)
0113     : DDDividedGeometryObject::DDDividedGeometryObject(div, cpv) {
0114   checkParametersValidity();
0115   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0116   setType("DivisionPolyconePhi");
0117   // this is the g4.  what do we keep??? I think it is deltaPhi
0118   // double deltaPhi = msol->GetEndPhi() - msol->GetStartPhi();
0119   if (divisionType_ == DivWIDTH) {
0120     //If you divide a tube of 360 degrees the offset displaces the starting angle, but you still fill the 360 degrees
0121     if (msol.deltaPhi() == 360._deg) {
0122       compNDiv_ = calculateNDiv(msol.deltaPhi(), div_.width(), 0.);
0123     } else {
0124       compNDiv_ = calculateNDiv(msol.deltaPhi(), div_.width(), div_.offset());
0125     }
0126   } else if (divisionType_ == DivNDIV) {
0127     if (msol.deltaPhi() == 360._deg) {
0128       compWidth_ = calculateWidth(msol.deltaPhi(), div_.nReplicas(), 0.);
0129     } else {
0130       compWidth_ = calculateWidth(msol.deltaPhi(), div_.nReplicas(), div_.offset());
0131     }
0132   }
0133 }
0134
0135 void DDDividedPolyconePhi::checkParametersValidity(void) { DDDividedGeometryObject::checkParametersValidity(); }
0136
0137 double DDDividedPolyconePhi::getMaxParameter(void) const {
0138   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0139   return msol.deltaPhi();
0140 }
0141
0142 DDRotation DDDividedPolyconePhi::makeDDRotation(const int copyNo) const {
0143   DDRotation myddrot;  // sets to identity.
0144   double posi = (copyNo - 1) * compWidth_;
0145   // how to name the rotation??
0146   // i do not like this
0147   DDName ddrotname(div_.parent().ddname().name() + "_DIVCHILD_ROT" + std::to_string(copyNo),
0148                    div_.parent().ddname().ns());
0149   myddrot = DDrot(ddrotname, changeRotMatrix(posi));
0150
0151   return myddrot;
0152 }
0153
0154 DDTranslation DDDividedPolyconePhi::makeDDTranslation(const int copyNo) const {
0155   DDTranslation translation;
0156   return translation;
0157 }
0158
0159 DDLogicalPart DDDividedPolyconePhi::makeDDLogicalPart(const int copyNo) const {
0160   DDSolid ddpolycone;
0161   DDMaterial usemat(div_.parent().material());
0162
0163   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0164   std::vector<double> localrMaxVec = msol.rMaxVec();
0165   std::vector<double> localrMinVec = msol.rMinVec();
0166   std::vector<double> localzVec = msol.zVec();
0167
0168   DDName solname(div_.parent().ddname().name() + "_DIVCHILD", div_.parent().ddname().ns());
0169   DDSolid sol(solname);
0170   if (!sol.isDefined().second) {
0171     ddpolycone = DDSolidFactory::polycone(
0172         solname, msol.startPhi() + div_.offset(), compWidth_, localzVec, localrMinVec, localrMaxVec);
0173   }
0174   DDLogicalPart ddlp(solname);
0175   if (!ddlp.isDefined().second) {
0176     ddlp = DDLogicalPart(solname, usemat, ddpolycone);
0177   }
0178
0179   return ddlp;
0180 }
0181
0182 DDDividedPolyconeZ::DDDividedPolyconeZ(const DDDivision& div, DDCompactView* cpv)
0183     : DDDividedGeometryObject::DDDividedGeometryObject(div, cpv) {
0184   checkParametersValidity();
0185   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0186   std::vector<double> localrMaxVec = msol.rMaxVec();
0187   std::vector<double> localrMinVec = msol.rMinVec();
0188   std::vector<double> localzVec = msol.zVec();
0189
0190   if (divisionType_ == DivWIDTH) {
0191     compNDiv_ = calculateNDiv(localzVec[localzVec.size() - 1] - localzVec[0], div_.width(), div_.offset());
0192   } else if (divisionType_ == DivNDIV) {
0193     compWidth_ = calculateNDiv(localzVec[localzVec.size() - 1] - localzVec[0], div_.nReplicas(), div_.offset());
0194   }
0195 }
0196
0197 void DDDividedPolyconeZ::checkParametersValidity(void) {
0198   DDDividedGeometryObject::checkParametersValidity();
0199
0200   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0201   std::vector<double> localzVec = msol.zVec();
0202   // CHANGE FROM G4 a polycone can be divided in Z by specifying
0203   // nReplicas IF they happen to coincide with the number of
0204   // z plans.
0205   size_t tempNDiv = div_.nReplicas();
0206   if (tempNDiv == 0)
0207     tempNDiv = calculateNDiv(localzVec[localzVec.size() - 1] - localzVec[0], div_.width(), div_.offset());
0208   if ((msol.zVec().size() - 1) != tempNDiv) {
0209     std::string s = "ERROR - DDDividedPolyconeZ::checkParametersValidity()";
0210     s += "\n\tDivision along Z will be done splitting in the defined";
0211     s += "\n\tz_planes, i.e, the number of division would be :";
0212     s += "\n\t" + std::to_string(msol.zVec().size() - 1);
0213     s += "\n\tinstead of " + std::to_string(tempNDiv) + " !\n";
0214
0215     throw cms::Exception("DDException") << s;
0216   }
0217 }
0218
0219 double DDDividedPolyconeZ::getMaxParameter(void) const {
0220   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0221   std::vector<double> localzVec = msol.zVec();
0222
0223   return (localzVec[localzVec.size() - 1] - localzVec[0]);
0224 }
0225
0226 DDRotation DDDividedPolyconeZ::makeDDRotation(const int copyNo) const {
0227   DDRotation myddrot;  // sets to identity.
0228   return myddrot;
0229 }
0230
0231 DDTranslation DDDividedPolyconeZ::makeDDTranslation(const int copyNo) const {
0232   DDTranslation translation;
0233   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0234   std::vector<double> localzVec = msol.zVec();
0235   double posi = (localzVec[copyNo] + localzVec[copyNo + 1]) / 2;
0236   translation.SetZ(posi);
0237   return translation;
0238 }
0239
0240 DDLogicalPart DDDividedPolyconeZ::makeDDLogicalPart(const int copyNo) const {
0241   DDSolid ddpolycone;
0242   DDMaterial usemat(div_.parent().material());
0243
0244   DDPolycone msol = (DDPolycone)(div_.parent().solid());
0245   std::vector<double> localrMaxVec = msol.rMaxVec();
0246   std::vector<double> localrMinVec = msol.rMinVec();
0247   std::vector<double> localzVec = msol.zVec();
0248
0249   DDName solname(div_.parent().ddname().name() + "_DIVCHILD" + std::to_string(copyNo), div_.parent().ddname().ns());
0250   ddpolycone = DDSolidFactory::cons(solname,
0251                                     compWidth_ / 2,
0252                                     localrMinVec[copyNo],
0253                                     localrMaxVec[copyNo],
0254                                     localrMinVec[copyNo + 1],
0255                                     localrMaxVec[copyNo + 1],
0256                                     msol.startPhi(),
0257                                     msol.deltaPhi());
0258
0259   DDLogicalPart ddlp = DDLogicalPart(solname, usemat, ddpolycone);
0260
0261   return ddlp;
0262 }