PrintGeomInfo/plugins/PrintG4Touch.cc

0001 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0002 #include "SimG4Core/Notification/interface/Observer.h"
0003 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0004 #include "SimG4Core/Geometry/interface/DD4hep2DDDName.h"
0005
0006 #include "FWCore/Framework/interface/EventSetup.h"
0007 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0008 #include "DataFormats/Math/interface/angle_units.h"
0009 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0010
0011 #include "G4LogicalVolumeStore.hh"
0012 #include "G4LogicalVolume.hh"
0013 #include "G4VSolid.hh"
0014 #include "G4Material.hh"
0015 #include "G4NavigationHistory.hh"
0016 #include "G4PhysicalVolumeStore.hh"
0017 #include "G4Region.hh"
0018 #include "G4RegionStore.hh"
0019 #include "G4Run.hh"
0020 #include "G4Track.hh"
0021 #include "G4TransportationManager.hh"
0022 #include "G4UserLimits.hh"
0023 #include "G4VisAttributes.hh"
0024 #include "G4VPhysicalVolume.hh"
0025
0026 #include <algorithm>
0027 #include <fstream>
0028 #include <map>
0029 #include <set>
0030 #include <string>
0031 #include <vector>
0032
0033 using angle_units::operators::convertRadToDeg;
0034
0035 class PrintG4Touch : public SimWatcher, public Observer<const BeginOfRun *> {
0036 public:
0037   PrintG4Touch(edm::ParameterSet const &p);
0038   ~PrintG4Touch() override = default;
0039
0040 private:
0041   void update(const BeginOfRun *run) override;
0042   void dumpSummary(std::ostream &out = G4cout);
0043   int countNoTouchables();
0044   void add1touchable(G4LogicalVolume *lv, int &nTouch);
0045   void dumpTouch(G4VPhysicalVolume *pv, unsigned int leafDepth, std::ostream &out = G4cout);
0046   void getTouch(G4VPhysicalVolume *pv, unsigned int leafDepth, unsigned int copym, std::vector<std::string> &touches);
0047   G4VPhysicalVolume *getTopPV();
0048   G4LogicalVolume *getTopLV();
0049
0050 private:
0051   bool dd4hep_, verbosity_;
0052   G4VPhysicalVolume *theTopPV_;
0053   G4NavigationHistory fHistory_;
0054 };
0055
0056 PrintG4Touch::PrintG4Touch(const edm::ParameterSet &p) {
0057   dd4hep_ = p.getUntrackedParameter<bool>("dd4hep", false);
0058   verbosity_ = p.getUntrackedParameter<bool>("verbosity", false);
0059   G4cout << "PrintG4Touch:: initialised for dd4hep " << dd4hep_ << " with verbosity levels:" << verbosity_ << G4endl;
0060 }
0061
0062 void PrintG4Touch::update(const BeginOfRun *run) {
0063   //Now take action
0064   theTopPV_ = getTopPV();
0065
0066   dumpSummary(G4cout);
0067
0068   std::vector<std::string> touches;
0069   getTouch(theTopPV_, 0, 1, touches);
0070   std::sort(touches.begin(), touches.end());
0071   for (const auto &touch : touches)
0072     G4cout << touch << G4endl;
0073
0074   //---------- Dump LV and PV information
0075   if (verbosity_)
0076     dumpTouch(theTopPV_, 0, G4cout);
0077 }
0078
0079 void PrintG4Touch::dumpSummary(std::ostream &out) {
0080   //---------- Dump number of objects of each class
0081   out << " @@@@@@@@@@@@@@@@@@ Dumping G4 geometry objects Summary " << G4endl;
0082   if (theTopPV_ == nullptr) {
0083     out << " No volume created " << G4endl;
0084     return;
0085   }
0086   out << " @@@ Geometry built inside world volume: "
0087       << DD4hep2DDDName::namePV(static_cast<std::string>(theTopPV_->GetName()), dd4hep_) << G4endl;
0088   // Get number of solids (< # LV if several LV share a solid)
0089   const G4LogicalVolumeStore *lvs = G4LogicalVolumeStore::GetInstance();
0090   std::vector<G4LogicalVolume *>::const_iterator lvcite;
0091   std::set<G4VSolid *> theSolids;
0092   for (lvcite = lvs->begin(); lvcite != lvs->end(); lvcite++)
0093     theSolids.insert((*lvcite)->GetSolid());
0094   out << " Number of G4VSolid's: " << theSolids.size() << G4endl;
0095   out << " Number of G4LogicalVolume's: " << lvs->size() << G4endl;
0096   const G4PhysicalVolumeStore *pvs = G4PhysicalVolumeStore::GetInstance();
0097   out << " Number of G4VPhysicalVolume's: " << pvs->size() << G4endl;
0098   out << " Number of Touchable's: " << countNoTouchables() << G4endl;
0099   const G4MaterialTable *matTab = G4Material::GetMaterialTable();
0100   out << " Number of G4Material's: " << matTab->size() << G4endl;
0101   const G4RegionStore *regs = G4RegionStore::GetInstance();
0102   out << " Number of G4Region's: " << regs->size() << G4endl;
0103 }
0104
0105 int PrintG4Touch::countNoTouchables() {
0106   int nTouch = 0;
0107   G4LogicalVolume *lv = getTopLV();
0108   add1touchable(lv, nTouch);
0109   return nTouch;
0110 }
0111
0112 void PrintG4Touch::add1touchable(G4LogicalVolume *lv, int &nTouch) {
0113   int siz = lv->GetNoDaughters();
0114   for (int ii = 0; ii < siz; ii++)
0115     add1touchable(lv->GetDaughter(ii)->GetLogicalVolume(), ++nTouch);
0116 }
0117
0118 void PrintG4Touch::dumpTouch(G4VPhysicalVolume *pv, unsigned int leafDepth, std::ostream &out) {
0119   if (leafDepth == 0)
0120     fHistory_.SetFirstEntry(pv);
0121   else
0122     fHistory_.NewLevel(pv, kNormal, pv->GetCopyNo());
0123
0124   G4LogicalVolume *lv = pv->GetLogicalVolume();
0125
0126   G4ThreeVector globalpoint = fHistory_.GetTopTransform().Inverse().TransformPoint(G4ThreeVector(0, 0, 0));
0127   std::string mother = (pv->GetMotherLogical())
0128                            ? (DD4hep2DDDName::nameSolid(
0129                                  static_cast<std::string>(pv->GetMotherLogical()->GetSolid()->GetName()), dd4hep_))
0130                            : "World";
0131   std::string lvname = DD4hep2DDDName::nameSolid(static_cast<std::string>(lv->GetName()), dd4hep_);
0132   out << leafDepth << "### VOLUME = " << lvname << " Copy No " << pv->GetCopyNo() << " in " << mother
0133       << " global position of centre " << globalpoint << " (r = " << globalpoint.perp()
0134       << ", phi = " << convertRadToDeg(globalpoint.phi()) << ")" << G4endl;
0135
0136   int NoDaughters = lv->GetNoDaughters();
0137   while ((NoDaughters--) > 0) {
0138     G4VPhysicalVolume *pvD = lv->GetDaughter(NoDaughters);
0139     if (!pvD->IsReplicated())
0140       dumpTouch(pvD, leafDepth + 1, out);
0141   }
0142
0143   if (leafDepth > 0)
0144     fHistory_.BackLevel();
0145 }
0146
0147 void PrintG4Touch::getTouch(G4VPhysicalVolume *pv,
0148                             unsigned int leafDepth,
0149                             unsigned int copym,
0150                             std::vector<std::string> &touches) {
0151   if (leafDepth == 0)
0152     fHistory_.SetFirstEntry(pv);
0153   else
0154     fHistory_.NewLevel(pv, kNormal, pv->GetCopyNo());
0155
0156   std::string mother = (pv->GetMotherLogical())
0157                            ? (DD4hep2DDDName::nameSolid(
0158                                  static_cast<std::string>(pv->GetMotherLogical()->GetSolid()->GetName()), dd4hep_))
0159                            : "World";
0160
0161   G4LogicalVolume *lv = pv->GetLogicalVolume();
0162   std::string lvname = DD4hep2DDDName::nameSolid(static_cast<std::string>(lv->GetSolid()->GetName()), dd4hep_);
0163   unsigned int copy = static_cast<unsigned int>(pv->GetCopyNo());
0164
0165   std::string type = static_cast<std::string>(lv->GetSolid()->GetEntityType());
0166
0167   std::string name = lvname + " " + std::to_string(copy) + " " + mother + " " + std::to_string(copym) + " " + type;
0168   touches.emplace_back(name);
0169
0170   int NoDaughters = lv->GetNoDaughters();
0171   while ((NoDaughters--) > 0) {
0172     G4VPhysicalVolume *pvD = lv->GetDaughter(NoDaughters);
0173     if (!pvD->IsReplicated())
0174       getTouch(pvD, leafDepth + 1, copy, touches);
0175   }
0176
0177   if (leafDepth > 0)
0178     fHistory_.BackLevel();
0179 }
0180
0181 G4VPhysicalVolume *PrintG4Touch::getTopPV() {
0182   return G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume();
0183 }
0184
0185 G4LogicalVolume *PrintG4Touch::getTopLV() { return theTopPV_->GetLogicalVolume(); }
0186
0187 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0188 #include "FWCore/PluginManager/interface/ModuleDef.h"
0189
0190 DEFINE_SIMWATCHER(PrintG4Touch);