File indexing completed on 2024-04-06 12:30:29
0001 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0002 #include "SimG4Core/Notification/interface/Observer.h"
0003 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0004
0005 #include "FWCore/Framework/interface/EventSetup.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "DataFormats/Math/interface/angle_units.h"
0008
0009 #include "G4Box.hh"
0010 #include "G4Cons.hh"
0011 #include "G4ExtrudedSolid.hh"
0012 #include "G4Polycone.hh"
0013 #include "G4Polyhedra.hh"
0014 #include "G4Trap.hh"
0015 #include "G4Trd.hh"
0016 #include "G4Tubs.hh"
0017 #include "G4LogicalVolumeStore.hh"
0018 #include "G4LogicalVolume.hh"
0019 #include "G4VSolid.hh"
0020 #include "G4NavigationHistory.hh"
0021 #include "G4PhysicalVolumeStore.hh"
0022 #include "G4TransportationManager.hh"
0023 #include "G4VPhysicalVolume.hh"
0024
0025 #include <algorithm>
0026 #include <map>
0027 #include <set>
0028 #include <sstream>
0029 #include <string>
0030 #include <vector>
0031
0032 using angle_units::operators::convertRadToDeg;
0033
0034 class PrintG4Solids : public SimWatcher, public Observer<const BeginOfRun *> {
0035 public:
0036 PrintG4Solids(edm::ParameterSet const &p);
0037 ~PrintG4Solids() override = default;
0038
0039 private:
0040 void update(const BeginOfRun *run) override;
0041 void dumpSummary(std::ostream &out = G4cout);
0042 G4VPhysicalVolume *getTopPV();
0043 bool select(const std::string &name, const std::string &shape) const;
0044 std::string reducedName(const std::string &name);
0045
0046 private:
0047 const bool dd4hep_;
0048 const std::vector<std::string> solids_;
0049 const std::vector<std::string> types_;
0050 G4VPhysicalVolume *theTopPV_;
0051 };
0052
0053 PrintG4Solids::PrintG4Solids(const edm::ParameterSet &p)
0054 : dd4hep_(p.getUntrackedParameter<bool>("dd4hep")),
0055 solids_(p.getUntrackedParameter<std::vector<std::string> >("dumpVolumes")),
0056 types_(p.getUntrackedParameter<std::vector<std::string> >("dumpShapes")) {
0057 G4cout << "PrintG4Solids:: initialised for printing information about G4VSolids for version dd4heP:" << dd4hep_
0058 << G4endl;
0059 }
0060
0061 void PrintG4Solids::update(const BeginOfRun *run) {
0062
0063 theTopPV_ = getTopPV();
0064
0065 dumpSummary(G4cout);
0066 }
0067
0068 void PrintG4Solids::dumpSummary(std::ostream &out) {
0069
0070 out << " @@@@@@@@@@@@@@@@@@ Dumping G4 geometry objects Summary " << G4endl;
0071 if (theTopPV_ == nullptr) {
0072 out << " No volume created " << G4endl;
0073 return;
0074 }
0075 out << " @@@ Geometry built inside world volume: " << theTopPV_->GetName() << G4endl;
0076
0077 const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
0078 std::vector<G4LogicalVolume *>::const_iterator lvcite;
0079 std::set<G4VSolid *> theSolids;
0080 for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++) {
0081 G4VSolid *solid = (*lvcite)->GetSolid();
0082 std::string name = static_cast<std::string>(solid->GetName());
0083 if (dd4hep_)
0084 name = reducedName(name);
0085 std::string type = static_cast<std::string>(solid->GetEntityType());
0086 if (select(name, type))
0087 theSolids.insert(solid);
0088 }
0089 out << " Number of G4VSolid's: " << theSolids.size() << G4endl;
0090 std::set<G4VSolid *>::const_iterator solid;
0091 for (solid = theSolids.begin(); solid != theSolids.end(); solid++) {
0092 G4String type = (*solid)->GetEntityType();
0093 std::string name = static_cast<std::string>((*solid)->GetName());
0094 if (dd4hep_)
0095 name = reducedName(name);
0096 out << name << ":" << type << " Volume " << (*solid)->GetCubicVolume();
0097 if (type == "G4Box") {
0098 const G4Box *box = static_cast<const G4Box *>(*solid);
0099 out << " dx:dy:dz " << box->GetXHalfLength() << ":" << box->GetYHalfLength() << ":" << box->GetZHalfLength();
0100 } else if (type == "G4Tubs") {
0101 const G4Tubs *tube = static_cast<const G4Tubs *>(*solid);
0102 out << " rin:rout:dz:phistart:dphi " << tube->GetInnerRadius() << ":" << tube->GetOuterRadius() << ":"
0103 << tube->GetZHalfLength() << ":" << convertRadToDeg(tube->GetStartPhiAngle()) << ":"
0104 << convertRadToDeg(tube->GetDeltaPhiAngle());
0105 } else if (type == "G4Cons") {
0106 const G4Cons *cone = static_cast<const G4Cons *>(*solid);
0107 out << " rinminus:routminus:rinplus:routplus:dz:phistart:dphi " << cone->GetInnerRadiusMinusZ() << ":"
0108 << cone->GetOuterRadiusMinusZ() << ":" << cone->GetInnerRadiusPlusZ() << ":" << cone->GetOuterRadiusPlusZ()
0109 << ":" << cone->GetZHalfLength() << ":" << convertRadToDeg(cone->GetStartPhiAngle()) << ":"
0110 << convertRadToDeg(cone->GetDeltaPhiAngle());
0111 } else if (type == "G4Trap") {
0112 const G4Trap *trap = static_cast<const G4Trap *>(*solid);
0113 out << " zhalf:yl1:xl11:xl12:tana1:yl2:xl21:xl22:tana2 " << trap->GetZHalfLength() << ":"
0114 << trap->GetYHalfLength1() << ":" << trap->GetXHalfLength1() << ":" << trap->GetXHalfLength2() << ":"
0115 << trap->GetTanAlpha1() << ":" << trap->GetYHalfLength2() << ":" << trap->GetXHalfLength3() << ":"
0116 << trap->GetXHalfLength4() << ":" << trap->GetTanAlpha2();
0117 } else if (type == "G4Trd") {
0118 const G4Trd *trd = static_cast<const G4Trd *>(*solid);
0119 out << " xl1:xl2:yl1:yl2:zhalf " << trd->GetXHalfLength1() << ":" << trd->GetXHalfLength2() << ":"
0120 << trd->GetYHalfLength1() << ":" << trd->GetYHalfLength2() << ":" << trd->GetZHalfLength();
0121 } else if (type == "G4Polycone") {
0122 const G4Polycone *cone = static_cast<const G4Polycone *>(*solid);
0123 const auto hist = cone->GetOriginalParameters();
0124 int num = hist->Num_z_planes;
0125 out << " angle " << convertRadToDeg(hist->Start_angle) << ":" << convertRadToDeg(hist->Opening_angle) << " with "
0126 << num << " planes:";
0127 for (int k = 0; k < num; ++k)
0128 out << " [" << k << "] " << hist->Z_values[k] << ":" << hist->Rmin[k] << ":" << hist->Rmax[k];
0129 } else if (type == "G4Polyhedra") {
0130 const G4Polyhedra *pgon = static_cast<const G4Polyhedra *>(*solid);
0131 const auto hist = pgon->GetOriginalParameters();
0132 int num = hist->Num_z_planes;
0133 out << " angle " << convertRadToDeg(hist->Start_angle) << ":" << convertRadToDeg(hist->Opening_angle) << " with "
0134 << hist->numSide << " sides and " << num << " planes:";
0135 for (int k = 0; k < num; ++k)
0136 out << " [" << k << "] " << hist->Z_values[k] << ":" << hist->Rmin[k] << ":" << hist->Rmax[k];
0137 } else if (type == "G4ExtrudedSolid") {
0138 const G4ExtrudedSolid *pgon = static_cast<const G4ExtrudedSolid *>(*solid);
0139 int vert = pgon->GetNofVertices();
0140 int numz = pgon->GetNofZSections();
0141 out << " " << vert << " vertices:";
0142 for (int k = 0; k < vert; ++k)
0143 out << " [" << k << "] " << pgon->GetVertex(k);
0144 out << "; and " << numz << " z-sections:";
0145 for (int k = 0; k < numz; ++k) {
0146 const auto &zsec = pgon->GetZSection(k);
0147 out << " [" << k << "] " << zsec.fZ << ":" << zsec.fScale << ":" << zsec.fOffset;
0148 }
0149 }
0150 out << G4endl;
0151 }
0152 }
0153
0154 G4VPhysicalVolume *PrintG4Solids::getTopPV() {
0155 return G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume();
0156 }
0157
0158 std::string PrintG4Solids::reducedName(const std::string &name) {
0159 std::string nam(name);
0160 uint32_t first = ((name.find(':') == std::string::npos) ? 0 : (name.find(':') + 1));
0161 uint32_t last(name.size() + 1);
0162 uint32_t loc(first);
0163 while (true) {
0164 if (name.find('_', loc) == std::string::npos)
0165 break;
0166 if (((loc + 5) < name.size()) && (name.substr(loc, 5) == "shape")) {
0167 last = loc;
0168 break;
0169 }
0170 loc = name.find('_', loc) + 1;
0171 if (loc > name.size())
0172 break;
0173 }
0174 nam = name.substr(first, last - first - 1);
0175 if ((last < name.size()) && (name.substr(name.size() - 5, 5) == "_refl"))
0176 nam += "_refl";
0177 return nam;
0178 }
0179
0180 bool PrintG4Solids::select(const std::string &name, const std::string &type) const {
0181 bool flag(true);
0182 if (!solids_.empty())
0183 flag = (flag && (std::find(solids_.begin(), solids_.end(), name) != solids_.end()));
0184 if (!types_.empty())
0185 flag = (flag && (std::find(types_.begin(), types_.end(), type) != types_.end()));
0186 return flag;
0187 }
0188
0189 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0190 #include "FWCore/PluginManager/interface/ModuleDef.h"
0191
0192 DEFINE_SIMWATCHER(PrintG4Solids);