Geometry/src/DDG4Builder.cc

0001 #include "FWCore/Utilities/interface/Exception.h"
0002
0003 #include "SimG4Core/Geometry/interface/DDG4Builder.h"
0004 #include "SimG4Core/Geometry/interface/DDG4SensitiveConverter.h"
0005 #include "SimG4Core/Geometry/interface/DDG4SolidConverter.h"
0006 #include "SimG4Core/Geometry/interface/SensitiveDetectorCatalog.h"
0007
0008 #include "DetectorDescription/Core/interface/DDCompactView.h"
0009 #include "DetectorDescription/Core/interface/DDSpecifics.h"
0010
0011 #include "G4LogicalVolume.hh"
0012 #include "G4Material.hh"
0013 #include "G4PVPlacement.hh"
0014 #include "G4ReflectionFactory.hh"
0015 #include "G4VPhysicalVolume.hh"
0016 #include "G4VSolid.hh"
0017
0018 #include <CLHEP/Units/SystemOfUnits.h>
0019 #include "G4UnitsTable.hh"
0020
0021 #include <sstream>
0022
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024
0025 DDG4Builder::DDG4Builder(const DDCompactView *cpv, G4LogicalVolumeToDDLogicalPartMap &lvmap, bool check)
0026     : solidConverter_(new DDG4SolidConverter), compactView_(cpv), map_(lvmap), check_(check) {
0027   theVectorOfDDG4Dispatchables_ = new DDG4DispContainer();
0028 }
0029
0030 DDG4Builder::~DDG4Builder() { delete solidConverter_; }
0031
0032 G4LogicalVolume *DDG4Builder::convertLV(const DDLogicalPart &part) {
0033   edm::LogVerbatim("SimG4CoreGeometry") << "DDG4Builder::convertLV(): DDLogicalPart = " << part;
0034   G4LogicalVolume *result = logs_[part];
0035   if (!result) {
0036     G4VSolid *g4s = convertSolid(part.solid());
0037     G4Material *g4m = convertMaterial(part.material());
0038     result = new G4LogicalVolume(g4s, g4m, part.name().name());
0039     map_.insert(result, part);
0040     DDG4Dispatchable *disp = new DDG4Dispatchable(&part, result);
0041     theVectorOfDDG4Dispatchables_->push_back(disp);
0042     edm::LogVerbatim("SimG4CoreGeometry")
0043         << "DDG4Builder::convertLV(): new G4LogicalVolume " << part.name().name()
0044         << "\nDDG4Builder: newEvent: dd=" << part.ddname() << " g4=" << result->GetName();
0045     logs_[part] = result;  // DDD -> GEANT4
0046   }
0047   return result;
0048 }
0049
0050 G4VSolid *DDG4Builder::convertSolid(const DDSolid &solid) {
0051   G4VSolid *result = sols_[solid];
0052   if (!result) {
0053     result = solidConverter_->convert(solid);
0054     sols_[solid] = result;
0055   }
0056   return result;
0057 }
0058
0059 G4Material *DDG4Builder::convertMaterial(const DDMaterial &material) {
0060   edm::LogVerbatim("SimG4CoreGeometry") << "DDDetConstr::ConvertMaterial: material=" << material;
0061   G4Material *result = nullptr;
0062   if (material) {
0063     // only if it's a valid DDD-material
0064     if ((result = mats_[material])) {
0065       edm::LogVerbatim("SimG4CoreGeometry") << "  is already converted";
0066       return result;
0067     }
0068   } else {
0069     // only if it's NOT a valid DDD-material
0070     throw cms::Exception("SimG4CoreGeometry",
0071                          " material is not valid from the Detector Description: " + material.toString());
0072   }
0073   int c = 0;
0074   if ((c = material.noOfConstituents())) {
0075     // it's a composite material
0076     edm::LogVerbatim("SimG4CoreGeometry")
0077         << "  creating a G4-composite material. c=" << c << " d=" << material.density() / CLHEP::g * CLHEP::mole;
0078     result = new G4Material(material.name().name(), material.density(), c);
0079     for (int i = 0; i < c; ++i) {
0080       // recursive building of constituents
0081       edm::LogVerbatim("SimG4CoreGeometry")
0082           << "  adding the composite=" << material.name() << " fm=" << material.constituent(i).second;
0083       result->AddMaterial(convertMaterial(material.constituent(i).first),
0084                           material.constituent(i).second);  // fractionmass
0085     }
0086   } else {
0087     // it's an elementary material
0088     edm::LogVerbatim("SimG4CoreGeometry") << "  building an elementary material"
0089                                           << " z=" << material.z() << " a=" << material.a() / CLHEP::g * CLHEP::mole
0090                                           << " d=" << material.density() / CLHEP::g * CLHEP::cm3;
0091     result = new G4Material(material.name().name(), material.z(), material.a(), material.density());
0092   }
0093   mats_[material] = result;
0094   return result;
0095 }
0096
0097 G4LogicalVolume *DDG4Builder::BuildGeometry(SensitiveDetectorCatalog &catalog) {
0098   G4ReflectionFactory *refFact = G4ReflectionFactory::Instance();
0099   refFact->SetScalePrecision(100. * refFact->GetScalePrecision());
0100
0101   using Graph = DDCompactView::Graph;
0102   const auto &gra = compactView_->graph();
0103   using adjl_iterator = Graph::const_adj_iterator;
0104   adjl_iterator git = gra.begin();
0105   adjl_iterator gend = gra.end();
0106
0107   for (; git != gend; ++git) {
0108     const DDLogicalPart &ddLP = gra.nodeData(git);
0109     if (!(ddLP.isDefined().second)) {
0110       throw cms::Exception("SimG4CoreGeometry",
0111                            " DDG4Builder::BuildGeometry() has encountered an "
0112                            "undefined DDLogicalPart named " +
0113                                ddLP.toString());
0114     }
0115     G4LogicalVolume *g4LV = convertLV(ddLP);
0116     if (!git->empty()) {
0117       // ask for children of ddLP
0118       Graph::edge_list::const_iterator cit = git->begin();
0119       Graph::edge_list::const_iterator cend = git->end();
0120       for (; cit != cend; ++cit) {
0121         // fetch specific data
0122         const DDLogicalPart &ddcurLP = gra.nodeData(cit->first);
0123         if (!ddcurLP.isDefined().second) {
0124           std::string err = " DDG4Builder::BuildGeometry() in processing \"children\" has ";
0125           err += "encountered an undefined DDLogicalPart named " + ddcurLP.toString() + " is a child of " +
0126                  ddLP.toString();
0127           throw cms::Exception("SimG4CoreGeometry", err);
0128         }
0129         int offset = getInt("CopyNoOffset", ddcurLP);
0130         int tag = getInt("CopyNoTag", ddcurLP);
0131         DDRotationMatrix rm(gra.edgeData(cit->second)->rot());
0132         DD3Vector x, y, z;
0133         rm.GetComponents(x, y, z);
0134         if ((x.Cross(y)).Dot(z) < 0)
0135           edm::LogVerbatim("SimG4CoreGeometry")
0136               << "DDG4Builder: Reflection: " << gra.edgeData(cit->second)->ddrot()
0137               << ">>Placement d=" << gra.nodeData(cit->first).ddname() << " m=" << ddLP.ddname()
0138               << " cp=" << gra.edgeData(cit->second)->copyno() << " r=" << gra.edgeData(cit->second)->ddrot().ddname();
0139         G4ThreeVector tempTran(gra.edgeData(cit->second)->trans().X(),
0140                                gra.edgeData(cit->second)->trans().Y(),
0141                                gra.edgeData(cit->second)->trans().Z());
0142         G4Translate3D transl = tempTran;
0143         CLHEP::HepRep3x3 temp(x.X(), x.Y(), x.Z(), y.X(), y.Y(), y.Z(), z.X(), z.Y(), z.Z());  // matrix
0144         CLHEP::HepRotation hr(temp);
0145         edm::LogVerbatim("SimG4CoreGeometry")
0146             << "Position " << gra.nodeData(cit->first).name().name() << ":"
0147             << gra.edgeData(cit->second)->copyno() + offset + tag << " in " << g4LV->GetName() << " at " << tempTran
0148             << " with rotation matrix (" << x.X() << ", " << x.Y() << ", " << x.Z() << ", " << y.X() << ", " << y.Y()
0149             << ", " << y.Z() << ", " << z.X() << ", " << z.Y() << ", " << z.Z() << ")";
0150
0151         // G3 convention of defining rot-matrices ...
0152         G4Transform3D trfrm = transl * G4Rotate3D(hr.inverse());  //.inverse();
0153
0154         refFact->Place(trfrm,  // transformation containing a possible reflection
0155                        gra.nodeData(cit->first).name().name(),
0156                        convertLV(gra.nodeData(cit->first)),                 // daugther
0157                        g4LV,                                                // mother
0158                        false,                                               // 'ONLY'
0159                        gra.edgeData(cit->second)->copyno() + offset + tag,  // copy number
0160                        check_);
0161       }  // iterate over children
0162     }    // if (children)
0163   }      // iterate over graph nodes
0164
0165   // Looking for in the G4ReflectionFactory secretly created reflected
0166   // G4LogicalVolumes
0167   std::map<DDLogicalPart, G4LogicalVolume *>::const_iterator ddg4_it = logs_.begin();
0168   for (; ddg4_it != logs_.end(); ++ddg4_it) {
0169     G4LogicalVolume *reflLogicalVolume = refFact->GetReflectedLV(ddg4_it->second);
0170     if (reflLogicalVolume) {
0171       DDLogicalPart ddlv = ddg4_it->first;
0172       map_.insert(reflLogicalVolume, ddlv);
0173       DDG4Dispatchable *disp = new DDG4Dispatchable(&(ddg4_it->first), reflLogicalVolume);
0174       theVectorOfDDG4Dispatchables_->push_back(disp);
0175       edm::LogVerbatim("SimG4CoreGeometry")
0176           << "DDG4Builder: dd=" << ddlv.ddname() << " g4=" << reflLogicalVolume->GetName();
0177     }
0178   }
0179
0180   G4LogicalVolume *world = logs_[compactView_->root()];
0181
0182   //
0183   //  needed for building sensitive detectors
0184   //
0185   DDG4SensitiveConverter conv;
0186   conv.upDate(*theVectorOfDDG4Dispatchables_, catalog);
0187
0188   return world;
0189 }
0190
0191 int DDG4Builder::getInt(const std::string &ss, const DDLogicalPart &part) {
0192   DDValue val(ss);
0193   std::vector<const DDsvalues_type *> result = part.specifics();
0194   bool foundIt = false;
0195   for (auto stype : result) {
0196     foundIt = DDfetch(stype, val);
0197     if (foundIt)
0198       break;
0199   }
0200   if (foundIt) {
0201     std::vector<double> temp = val.doubles();
0202     if (temp.size() != 1) {
0203       throw cms::Exception("SimG4CoreGeometry",
0204                            " DDG4Builder::getInt() Problem with Region tags - "
0205                            "one and only one allowed: " +
0206                                ss);
0207     }
0208     return int(temp[0]);
0209   } else
0210     return 0;
0211 }
0212
0213 double DDG4Builder::getDouble(const std::string &ss, const DDLogicalPart &part) {
0214   DDValue val(ss);
0215   std::vector<const DDsvalues_type *> result = part.specifics();
0216   bool foundIt = false;
0217   for (auto stype : result) {
0218     foundIt = DDfetch(stype, val);
0219     if (foundIt)
0220       break;
0221   }
0222   if (foundIt) {
0223     std::vector<std::string> temp = val.strings();
0224     if (temp.size() != 1) {
0225       throw cms::Exception("SimG4CoreGeometry",
0226                            " DDG4Builder::getDouble() Problem with Region tags "
0227                            "- one and only one allowed: " +
0228                                ss);
0229     }
0230     double v;
0231     std::string unit;
0232     std::istringstream is(temp[0]);
0233     is >> v >> unit;
0234     v = v * G4UnitDefinition::GetValueOf(unit.substr(1, unit.size()));
0235     return v;
0236   } else
0237     return 0;
0238 }