Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-05-10 02:21:25 

0001 #include "G4VSolid.hh"
0002 #include "SimG4Core/Geometry/interface/DDG4SolidConverter.h"
0003 
0004 #include "DetectorDescription/Core/interface/DDSolid.h"
0005 
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007 #include "FWCore/Utilities/interface/Exception.h"
0008 #include "G4ios.hh"
0009 #include <CLHEP/Units/SystemOfUnits.h>
0010 
0011 //#define EDM_ML_DEBUG
0012 
0013 const std::vector<double> *DDG4SolidConverter::par_ = nullptr;
0014 G4RotationMatrix *DDG4SolidConverter::rot = nullptr;
0015 
0016 DDG4SolidConverter::DDG4SolidConverter() {
0017   // could also be done 'dynamically' from outside
0018   // would then need to have a 'register' method ...
0019   convDispatch_[DDSolidShape::ddbox] = DDG4SolidConverter::box;
0020   convDispatch_[DDSolidShape::ddtubs] = DDG4SolidConverter::tubs;
0021   convDispatch_[DDSolidShape::ddcuttubs] = DDG4SolidConverter::cuttubs;
0022   convDispatch_[DDSolidShape::ddtrap] = DDG4SolidConverter::trap;
0023   convDispatch_[DDSolidShape::ddcons] = DDG4SolidConverter::cons;
0024   convDispatch_[DDSolidShape::ddpolycone_rrz] = DDG4SolidConverter::polycone_rrz;
0025   convDispatch_[DDSolidShape::ddpolycone_rz] = DDG4SolidConverter::polycone_rz;
0026   convDispatch_[DDSolidShape::ddpolyhedra_rrz] = DDG4SolidConverter::polyhedra_rrz;
0027   convDispatch_[DDSolidShape::ddpolyhedra_rz] = DDG4SolidConverter::polyhedra_rz;
0028   convDispatch_[DDSolidShape::ddextrudedpolygon] = DDG4SolidConverter::extrudedpolygon;
0029   convDispatch_[DDSolidShape::ddtorus] = DDG4SolidConverter::torus;
0030   convDispatch_[DDSolidShape::ddunion] = DDG4SolidConverter::unionsolid;
0031   convDispatch_[DDSolidShape::ddintersection] = DDG4SolidConverter::intersection;
0032   convDispatch_[DDSolidShape::ddsubtraction] = DDG4SolidConverter::subtraction;
0033   convDispatch_[DDSolidShape::ddpseudotrap] = DDG4SolidConverter::pseudotrap;
0034   convDispatch_[DDSolidShape::ddtrunctubs] = DDG4SolidConverter::trunctubs;
0035   convDispatch_[DDSolidShape::ddsphere] = DDG4SolidConverter::sphere;
0036   convDispatch_[DDSolidShape::ddellipticaltube] = DDG4SolidConverter::ellipticaltube;
0037 }
0038 
0039 DDG4SolidConverter::~DDG4SolidConverter() {}
0040 
0041 G4VSolid *DDG4SolidConverter::convert(const DDSolid &solid) {
0042   if (!solid) {
0043     edm::LogError("SimG4CoreGeometry") << " DDG4SolidConverter::convert(..) found an undefined DDSolid "
0044                                        << solid.toString();
0045     throw cms::Exception("SimG4CoreGeometry",
0046                          "DDG4SolidConverter::convert(..) found an undefined DDSolid " + solid.toString());
0047   }
0048   G4VSolid *result = nullptr;
0049   par_ = &(solid.parameters());
0050   std::map<DDSolidShape, FNPTR>::iterator it = convDispatch_.find(solid.shape());
0051   if (it != convDispatch_.end()) {
0052     result = it->second(solid);
0053   } else {
0054     throw cms::Exception("DetectorDescriptionFault")
0055         << "DDG4SolidConverter::convert: conversion failed for s=" << solid
0056         << "\n solid.shape()=" << DDSolidShapesName::name(solid.shape()) << G4endl;
0057   }
0058   return result;
0059 }
0060 
0061 #include "G4Box.hh"
0062 G4VSolid *DDG4SolidConverter::box(const DDSolid &solid) {
0063   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: box = " << solid;
0064   return new G4Box(solid.name().name(), (*par_)[0], (*par_)[1], (*par_)[2]);
0065 }
0066 
0067 #include "G4Tubs.hh"
0068 G4VSolid *DDG4SolidConverter::tubs(const DDSolid &solid) {
0069   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: tubs = " << solid;
0070   return new G4Tubs(solid.name().name(),
0071                     (*par_)[1],   // rmin
0072                     (*par_)[2],   // rmax
0073                     (*par_)[0],   // dzHalf
0074                     (*par_)[3],   // phiStart
0075                     (*par_)[4]);  // deltaPhi
0076 }
0077 
0078 #include "G4CutTubs.hh"
0079 G4VSolid *DDG4SolidConverter::cuttubs(const DDSolid &solid) {
0080   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: tubs = " << solid;
0081   return new G4CutTubs(solid.name().name(),
0082                        (*par_)[1],  // rmin
0083                        (*par_)[2],  // rmax
0084                        (*par_)[0],  // dzHalf
0085                        (*par_)[3],  // phiStart
0086                        (*par_)[4],  // deltaPhi
0087                        G4ThreeVector((*par_)[5], (*par_)[6], (*par_)[7]),
0088                        G4ThreeVector((*par_)[8], (*par_)[9], (*par_)[10]));
0089 }
0090 
0091 #include "G4Trap.hh"
0092 G4VSolid *DDG4SolidConverter::trap(const DDSolid &solid) {
0093   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: trap = " << solid;
0094   return new G4Trap(solid.name().name(),
0095                     (*par_)[0],    // pDz
0096                     (*par_)[1],    // theta
0097                     (*par_)[2],    // phi
0098                     (*par_)[3],    // y1
0099                     (*par_)[4],    // x1
0100                     (*par_)[5],    // x2
0101                     (*par_)[6],    // alpha1
0102                     (*par_)[7],    // y2
0103                     (*par_)[8],    // x3
0104                     (*par_)[9],    // x4
0105                     (*par_)[10]);  // alpha2
0106 }
0107 
0108 #include "G4Cons.hh"
0109 G4VSolid *DDG4SolidConverter::cons(const DDSolid &solid) {
0110   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: cons = " << solid;
0111   return new G4Cons(solid.name().name(),
0112                     (*par_)[1],   // rmin -z
0113                     (*par_)[2],   // rmax -z
0114                     (*par_)[3],   // rmin +z
0115                     (*par_)[4],   // rmax +z
0116                     (*par_)[0],   // zHalf
0117                     (*par_)[5],   // phistart
0118                     (*par_)[6]);  // deltaphi
0119 }
0120 
0121 #include "G4Polycone.hh"
0122 G4VSolid *DDG4SolidConverter::polycone_rz(const DDSolid &solid) {
0123   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: pcon_rz = " << solid;
0124   std::vector<double> r;
0125   std::vector<double> z;
0126   std::vector<double>::const_iterator i = (*par_).begin() + 2;
0127   int count = 0;
0128   for (; i != (*par_).end(); ++i) {
0129     edm::LogVerbatim("SimG4CoreGeometry") << " z=" << *i / CLHEP::cm;
0130     z.push_back(*i);
0131     ++i;
0132     edm::LogVerbatim("SimG4CoreGeometry") << " r=" << *i / CLHEP::cm;
0133     r.push_back(*i);
0134     count++;
0135   }
0136   edm::LogVerbatim("SimG4CoreGeometry") << "sp=" << (*par_)[0] / CLHEP::deg << " ep=" << (*par_)[1] / CLHEP::deg;
0137 #ifdef EDM_ML_DEBUG
0138   G4cout << "### Polycone_RZ: "
0139          << "sp=" << (*par_)[0] / CLHEP::deg << " ep=" << (*par_)[1] / CLHEP::deg << " N= " << count << G4endl;
0140   for (int i = 0; i < count; ++i) {
0141     G4cout << " R= " << r[i] << " Z= " << z[i] << G4endl;
0142   }
0143 #endif
0144   return new G4Polycone(solid.name().name(),
0145                         (*par_)[0],
0146                         (*par_)[1],  // start,delta-phi
0147                         count,       // numRZ
0148                         &(r[0]),
0149                         &(z[0]));
0150 }
0151 
0152 G4VSolid *DDG4SolidConverter::polycone_rrz(const DDSolid &solid) {
0153   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: pcon_rrz = " << solid;
0154   std::vector<double> z_p;
0155   std::vector<double> rmin_p;
0156   std::vector<double> rmax_p;
0157   std::vector<double>::const_iterator i = par_->begin() + 2;
0158   int count = 0;
0159   for (; i != par_->end(); ++i) {
0160     edm::LogVerbatim("SimG4CoreGeometry") << "z=" << *i / CLHEP::cm;
0161     z_p.push_back(*i);
0162     ++i;
0163     edm::LogVerbatim("SimG4CoreGeometry") << "rmin=" << *i / CLHEP::cm;
0164     rmin_p.push_back(*i);
0165     ++i;
0166     edm::LogVerbatim("SimG4CoreGeometry") << "rmax=" << *i / CLHEP::cm;
0167     rmax_p.push_back(*i);
0168     count++;
0169   }
0170   edm::LogVerbatim("SimG4CoreGeometry") << "sp=" << (*par_)[0] / CLHEP::deg << " ep=" << (*par_)[1] / CLHEP::deg;
0171 #ifdef EDM_ML_DEBUG
0172   G4cout << "### Polycone_RRZ: "
0173          << "sp=" << (*par_)[0] / CLHEP::deg << " ep=" << (*par_)[1] / CLHEP::deg << " N= " << count << G4endl;
0174   for (int i = 0; i < count; ++i) {
0175     G4cout << " R1= " << rmin_p[i] << " R1= " << rmax_p[i] << " Z= " << z_p[i] << G4endl;
0176   }
0177 #endif
0178   return new G4Polycone(solid.name().name(),
0179                         (*par_)[0],
0180                         (*par_)[1],  // start,delta-phi
0181                         count,       // sections
0182                         &(z_p[0]),
0183                         &(rmin_p[0]),
0184                         &(rmax_p[0]));
0185 }
0186 
0187 #include "G4Polyhedra.hh"
0188 G4VSolid *DDG4SolidConverter::polyhedra_rz(const DDSolid &solid) {
0189   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: phed_rz = " << solid;
0190   std::vector<double> r;
0191   std::vector<double> z;
0192   std::vector<double>::const_iterator i = par_->begin() + 3;
0193   int count = 0;
0194 
0195   for (; i != par_->end(); ++i) {
0196     z.push_back(*i);
0197     ++i;
0198     r.push_back(*i);
0199     count++;
0200   }
0201 
0202   return new G4Polyhedra(solid.name().name(),
0203                          (*par_)[1],
0204                          (*par_)[2],
0205                          int((*par_)[0]),  // start,delta-phi;sides
0206                          count,            // numRZ
0207                          &(r[0]),
0208                          &(z[0]));
0209 }
0210 
0211 G4VSolid *DDG4SolidConverter::polyhedra_rrz(const DDSolid &solid) {
0212   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: phed_rrz = " << solid;
0213   std::vector<double> z_p;
0214   std::vector<double> rmin_p;
0215   std::vector<double> rmax_p;
0216   std::vector<double>::const_iterator i = par_->begin() + 3;
0217   int count = 0;
0218   for (; i != par_->end(); ++i) {
0219     edm::LogVerbatim("SimG4CoreGeometry") << "z=" << *i / CLHEP::cm;
0220     z_p.push_back(*i);
0221     ++i;
0222     edm::LogVerbatim("SimG4CoreGeometry") << "rmin=" << *i / CLHEP::cm;
0223     rmin_p.push_back(*i);
0224     ++i;
0225     edm::LogVerbatim("SimG4CoreGeometry") << "rmax=" << *i / CLHEP::cm;
0226     rmax_p.push_back(*i);
0227     count++;
0228   }
0229   edm::LogVerbatim("SimG4CoreGeometry") << "sp=" << (*par_)[0] / CLHEP::deg << " ep=" << (*par_)[1] / CLHEP::deg;
0230   return new G4Polyhedra(solid.name().name(),
0231                          (*par_)[1],
0232                          (*par_)[2],
0233                          int((*par_)[0]),  // start,delta-phi,sides
0234                          count,            // sections
0235                          &(z_p[0]),
0236                          &(rmin_p[0]),
0237                          &(rmax_p[0]));
0238 }
0239 
0240 #include "G4ExtrudedSolid.hh"
0241 G4VSolid *DDG4SolidConverter::extrudedpolygon(const DDSolid &solid) {
0242   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: extr_pgon = " << solid;
0243   std::vector<double> x = static_cast<DDExtrudedPolygon>(solid).xVec();
0244   std::vector<double> y = static_cast<DDExtrudedPolygon>(solid).yVec();
0245   std::vector<double> z = static_cast<DDExtrudedPolygon>(solid).zVec();
0246   std::vector<double> zx = static_cast<DDExtrudedPolygon>(solid).zxVec();
0247   std::vector<double> zy = static_cast<DDExtrudedPolygon>(solid).zyVec();
0248   std::vector<double> zs = static_cast<DDExtrudedPolygon>(solid).zscaleVec();
0249 
0250   std::vector<G4TwoVector> polygon;
0251   std::vector<G4ExtrudedSolid::ZSection> zsections;
0252   for (unsigned int it = 0; it < x.size(); ++it)
0253     polygon.emplace_back(x[it], y[it]);
0254   for (unsigned int it = 0; it < z.size(); ++it)
0255     zsections.emplace_back(z[it], G4TwoVector(zx[it], zy[it]), zs[it]);
0256   return new G4ExtrudedSolid(solid.name().name(), polygon, zsections);
0257 }
0258 
0259 #include "G4Torus.hh"
0260 G4VSolid *DDG4SolidConverter::torus(const DDSolid &solid) {
0261   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: torus = " << solid;
0262   return new G4Torus(solid.name().name(),
0263                      (*par_)[0],   // rmin
0264                      (*par_)[1],   // rmax
0265                      (*par_)[2],   // Rtor
0266                      (*par_)[3],   // phiStart
0267                      (*par_)[4]);  // deltaPhi
0268 }
0269 
0270 #include "G4UnionSolid.hh"
0271 G4VSolid *DDG4SolidConverter::unionsolid(const DDSolid &solid) {
0272   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: unionsolid = " << solid.name();
0273   G4UnionSolid *us = nullptr;
0274   DDBooleanSolid bs(solid);
0275   if (bs) {
0276     edm::LogVerbatim("SimG4CoreGeometry") << "SolidA=" << bs.solidA();
0277     G4VSolid *sa = DDG4SolidConverter().convert(bs.solidA());
0278     edm::LogVerbatim("SimG4CoreGeometry") << "SolidB=" << bs.solidB();
0279     G4VSolid *sb = DDG4SolidConverter().convert(bs.solidB());
0280     edm::LogVerbatim("SimG4CoreGeometry") << " name:" << solid.name() << " t=" << bs.translation() << std::flush;
0281     edm::LogVerbatim("SimG4CoreGeometry") << " " << bs.rotation().rotation().Inverse() << std::flush;
0282     std::vector<double> tdbl(9);
0283     bs.rotation().rotation().Inverse().GetComponents(tdbl.begin(), tdbl.end());
0284     CLHEP::HepRep3x3 temprep(tdbl[0], tdbl[1], tdbl[2], tdbl[3], tdbl[4], tdbl[5], tdbl[6], tdbl[7], tdbl[8]);
0285     CLHEP::Hep3Vector temphvec(bs.translation().X(), bs.translation().Y(), bs.translation().Z());
0286     us = new G4UnionSolid(solid.name().name(), sa, sb, new CLHEP::HepRotation(temprep), temphvec);
0287 
0288   }  // else?
0289   return us;
0290 }
0291 
0292 #include "G4SubtractionSolid.hh"
0293 #include <sstream>
0294 G4VSolid *DDG4SolidConverter::subtraction(const DDSolid &solid) {
0295   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: subtraction = " << solid;
0296   G4SubtractionSolid *us = nullptr;
0297   DDBooleanSolid bs(solid);
0298   if (bs) {
0299     G4VSolid *sa = DDG4SolidConverter().convert(bs.solidA());
0300     G4VSolid *sb = DDG4SolidConverter().convert(bs.solidB());
0301     edm::LogVerbatim("SimG4CoreGeometry") << " name:" << solid.name() << " t=" << bs.translation() << std::flush;
0302     edm::LogVerbatim("SimG4CoreGeometry") << " " << bs.rotation().rotation().Inverse() << std::flush;
0303     std::vector<double> tdbl(9);
0304     bs.rotation().rotation().Inverse().GetComponents(tdbl.begin(), tdbl.end());
0305     CLHEP::HepRep3x3 temprep(tdbl[0], tdbl[1], tdbl[2], tdbl[3], tdbl[4], tdbl[5], tdbl[6], tdbl[7], tdbl[8]);
0306     CLHEP::Hep3Vector temphvec(bs.translation().X(), bs.translation().Y(), bs.translation().Z());
0307     us = new G4SubtractionSolid(solid.name().name(), sa, sb, new CLHEP::HepRotation(temprep), temphvec);
0308   }
0309   return us;
0310 }
0311 
0312 #include "G4IntersectionSolid.hh"
0313 G4VSolid *DDG4SolidConverter::intersection(const DDSolid &solid) {
0314   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: intersection = " << solid;
0315   G4IntersectionSolid *us = nullptr;
0316   DDBooleanSolid bs(solid);
0317   if (bs) {
0318     G4VSolid *sa = DDG4SolidConverter().convert(bs.solidA());
0319     G4VSolid *sb = DDG4SolidConverter().convert(bs.solidB());
0320     edm::LogVerbatim("SimG4CoreGeometry") << " name:" << solid.name() << " t=" << bs.translation() << std::flush;
0321     edm::LogVerbatim("SimG4CoreGeometry") << " " << bs.rotation().rotation().Inverse() << std::flush;
0322     std::vector<double> tdbl(9);
0323     bs.rotation().rotation().Inverse().GetComponents(tdbl.begin(), tdbl.end());
0324     CLHEP::HepRep3x3 temprep(tdbl[0], tdbl[1], tdbl[2], tdbl[3], tdbl[4], tdbl[5], tdbl[6], tdbl[7], tdbl[8]);
0325     CLHEP::Hep3Vector temphvec(bs.translation().X(), bs.translation().Y(), bs.translation().Z());
0326     us = new G4IntersectionSolid(solid.name().name(), sa, sb, new CLHEP::HepRotation(temprep), temphvec);
0327   }
0328   return us;
0329 }
0330 
0331 #include "G4Trd.hh"
0332 G4VSolid *DDG4SolidConverter::pseudotrap(const DDSolid &solid) {
0333   if (nullptr == rot) {
0334     rot = new G4RotationMatrix;
0335     rot->rotateX(90. * CLHEP::deg);
0336   }
0337 
0338   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: pseudoTrap = " << solid;
0339   G4Trd *trap = nullptr;
0340   G4Tubs *tubs = nullptr;
0341   G4VSolid *result = nullptr;
0342   DDPseudoTrap pt(solid);  // pt...PseudoTrap
0343   double r = pt.radius();
0344   bool atMinusZ = pt.atMinusZ();
0345   double x = 0;
0346   double h = 0;
0347   bool intersec = false;  // union or intersection solid
0348   if (pt.atMinusZ()) {
0349     x = pt.x1();  // tubs radius
0350   } else {
0351     x = pt.x2();  // tubs radius
0352   }
0353   double openingAngle = 2. * asin(x / std::abs(r));
0354   // trap = new G4Trd(solid.name().name(),
0355   double displacement = 0;
0356   double startPhi = 0;
0357   /* calculate the displacement of the tubs w.r.t. to the trap,
0358      determine the opening angle of the tubs */
0359   double delta = sqrt(r * r - x * x);
0360   if (r < 0 && std::abs(r) >= x) {
0361     intersec = true;                            // intersection solid
0362     h = pt.y1() < pt.y2() ? pt.y2() : pt.y1();  // tubs half height
0363     h += h / 20.;                               // enlarge a bit - for subtraction solid
0364     if (atMinusZ) {
0365       displacement = -pt.halfZ() - delta;
0366       startPhi = 270. * CLHEP::deg - openingAngle / 2.;
0367     } else {
0368       displacement = pt.halfZ() + delta;
0369       startPhi = 90. * CLHEP::deg - openingAngle / 2.;
0370     }
0371   } else if (r > 0 && std::abs(r) >= x) {
0372     if (atMinusZ) {
0373       displacement = -pt.halfZ() + delta;
0374       startPhi = 90. * CLHEP::deg - openingAngle / 2.;
0375       h = pt.y1();
0376     } else {
0377       displacement = pt.halfZ() - delta;
0378       startPhi = 270. * CLHEP::deg - openingAngle / 2.;
0379       h = pt.y2();
0380     }
0381   } else {
0382     throw cms::Exception("DetectorDescriptionFault", "Check parameters of the PseudoTrap! name=" + pt.name().name());
0383   }
0384   G4ThreeVector displ(0., 0.,
0385                       displacement);  // displacement of the tubs w.r.t. trap
0386   edm::LogVerbatim("SimG4CoreGeometry") << "DDSolidConverter::pseudotrap(): displacement=" << displacement
0387                                         << " openingAngle=" << openingAngle / CLHEP::deg << " x=" << x << " h=" << h;
0388 
0389   // Now create two solids (trd & tubs), and a boolean solid out of them
0390   std::string name = pt.name().name();
0391   trap = new G4Trd(name, pt.x1(), pt.x2(), pt.y1(), pt.y2(), pt.halfZ());
0392   tubs = new G4Tubs(name,
0393                     0.,           // rMin
0394                     std::abs(r),  // rMax
0395                     h,            // half height
0396                     startPhi,     // start angle
0397                     openingAngle);
0398   if (intersec) {
0399     result = new G4SubtractionSolid(name, trap, tubs, rot, displ);
0400   } else {
0401     G4VSolid *tubicCap = new G4SubtractionSolid(
0402         name, tubs, new G4Box(name, 1.1 * x, sqrt(r * r - x * x), 1.1 * h), nullptr, G4ThreeVector());
0403     result = new G4UnionSolid(name, trap, tubicCap, rot, displ);
0404   }
0405   return result;
0406 }
0407 
0408 G4VSolid *DDG4SolidConverter::trunctubs(const DDSolid &solid) {
0409   // truncated tube-section: a boolean subtraction solid:
0410   //                         from a tube-section a box is subtracted according
0411   //                         to the given parameters
0412   edm::LogVerbatim("SimG4CoreGeometry") << "MantisConverter: solidshape=" << DDSolidShapesName::name(solid.shape())
0413                                         << " " << solid;
0414   edm::LogVerbatim("SimG4CoreGeometry") << "before";
0415   DDTruncTubs tt(solid);
0416   edm::LogVerbatim("SimG4CoreGeometry") << "after";
0417   double rIn(tt.rIn()), rOut(tt.rOut()), zHalf(tt.zHalf()), startPhi(tt.startPhi()), deltaPhi(tt.deltaPhi()),
0418       cutAtStart(tt.cutAtStart()), cutAtDelta(tt.cutAtDelta());
0419   bool cutInside(bool(tt.cutInside()));
0420   std::string name = tt.name().name();
0421 
0422   // check the parameters
0423   if (rIn <= 0 || rOut <= 0 || cutAtStart <= 0 || cutAtDelta <= 0) {
0424     throw cms::Exception(
0425         "DetectorDescriptionFault",
0426         "TruncTubs " + std::string(tt.name().fullname()) + ": 0 <= rIn,cutAtStart,rOut,cutAtDelta,rOut violated!");
0427   }
0428   if (rIn >= rOut) {
0429     throw cms::Exception("DetectorDescriptionFault",
0430                          "TruncTubs " + std::string(tt.name().fullname()) + ": rIn<rOut violated!");
0431   }
0432   if (startPhi != 0.) {
0433     throw cms::Exception("DetectorDescriptionFault",
0434                          "TruncTubs " + std::string(tt.name().fullname()) + ": startPhi != 0 not supported!");
0435   }
0436 
0437   startPhi = 0.;
0438   double r(cutAtStart), R(cutAtDelta);
0439   G4VSolid *result(nullptr);
0440   G4VSolid *tubs = new G4Tubs(name, rIn, rOut, zHalf, startPhi, deltaPhi);
0441   edm::LogVerbatim("SimG4CoreGeometry") << "G4Tubs: " << rIn / CLHEP::cm << ' ' << rOut / CLHEP::cm << ' '
0442                                         << zHalf / CLHEP::cm << ' ' << startPhi / CLHEP::deg << ' '
0443                                         << deltaPhi / CLHEP::deg;
0444   edm::LogVerbatim("SimG4CoreGeometry") << solid;
0445   // length & hight of the box
0446   double boxX(rOut * std::sqrt(2) * 1.1),
0447       boxY(rOut * std::sqrt(2) * 1.1);  // exaggerate dimensions - does not matter, it's subtracted!
0448 
0449   // width of the box > width of the tubs
0450   double boxZ(1.1 * zHalf);
0451 
0452   // angle of the box w.r.t. tubs
0453   double cath = r - R * cos(deltaPhi);
0454   double hypo = sqrt(r * r + R * R - 2. * r * R * cos(deltaPhi));
0455   double cos_alpha = cath / hypo;
0456 
0457   double alpha = -acos(cos_alpha);
0458   edm::LogVerbatim("SimG4CoreGeometry") << "cath=" << cath / CLHEP::cm;
0459   edm::LogVerbatim("SimG4CoreGeometry") << "hypo=" << hypo / CLHEP::cm;
0460   edm::LogVerbatim("SimG4CoreGeometry") << "al=" << acos(cath / hypo) / CLHEP::deg;
0461   edm::LogVerbatim("SimG4CoreGeometry") << "deltaPhi=" << deltaPhi / CLHEP::deg << "\n"
0462                                         << "r=" << r / CLHEP::cm << "\n"
0463                                         << "R=" << R / CLHEP::cm;
0464 
0465   edm::LogVerbatim("SimG4CoreGeometry") << "alpha=" << alpha / CLHEP::deg;
0466 
0467   // rotationmatrix of box w.r.t. tubs
0468   G4RotationMatrix *rot = new G4RotationMatrix;
0469   rot->rotateZ(-alpha);
0470   edm::LogVerbatim("SimG4CoreGeometry") << (*rot);
0471 
0472   // center point of the box
0473   double xBox;
0474   if (!cutInside) {
0475     xBox = r + boxY / sin(std::abs(alpha));
0476   } else {
0477     xBox = -(boxY / sin(std::abs(alpha)) - r);
0478   }
0479 
0480   G4ThreeVector trans(xBox, 0., 0.);
0481   edm::LogVerbatim("SimG4CoreGeometry") << "trans=" << trans;
0482 
0483   G4VSolid *box = new G4Box(name, boxX, boxY, boxZ);
0484   result = new G4SubtractionSolid(name, tubs, box, rot, trans);
0485 
0486   return result;
0487 }
0488 
0489 #include "G4Sphere.hh"
0490 G4VSolid *DDG4SolidConverter::sphere(const DDSolid &solid) {
0491   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: sphere = " << solid;
0492   DDSphere sp(solid);
0493   return new G4Sphere(solid.name().name(),
0494                       sp.innerRadius(),
0495                       sp.outerRadius(),
0496                       sp.startPhi(),
0497                       sp.deltaPhi(),
0498                       sp.startTheta(),
0499                       sp.deltaTheta());
0500 }
0501 
0502 #include "G4EllipticalTube.hh"
0503 G4VSolid *DDG4SolidConverter::ellipticaltube(const DDSolid &solid) {
0504   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4SolidConverter: ellipticaltube = " << solid;
0505   DDEllipticalTube sp(solid);
0506   return new G4EllipticalTube(solid.name().name(), sp.xSemiAxis(), sp.ySemiAxis(), sp.zHeight());
0507 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!