Forward/plugins/PrintMTDSens.cc

0001 //#define EDM_ML_DEBUG
0002
0003 #include "SimG4Core/Notification/interface/BeginOfRun.h"
0004 #include "SimG4Core/Notification/interface/Observer.h"
0005 #include "SimG4Core/Watcher/interface/SimWatcher.h"
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0008
0009 #include "Geometry/MTDCommonData/interface/BTLNumberingScheme.h"
0010 #include "Geometry/MTDCommonData/interface/ETLNumberingScheme.h"
0011 #include "DataFormats/ForwardDetId/interface/BTLDetId.h"
0012 #include "DataFormats/ForwardDetId/interface/ETLDetId.h"
0013 #include "DataFormats/Math/interface/Rounding.h"
0014
0015 #include "G4Run.hh"
0016 #include "G4VPhysicalVolume.hh"
0017 #include "G4LogicalVolume.hh"
0018 #include "G4NavigationHistory.hh"
0019 #include "G4TransportationManager.hh"
0020 #include "G4Box.hh"
0021
0022 #include <set>
0023 #include <map>
0024 #include <iostream>
0025 #include <string>
0026 #include <vector>
0027
0028 class PrintMTDSens : public SimWatcher, public Observer<const BeginOfRun *> {
0029 public:
0030   PrintMTDSens(edm::ParameterSet const &p);
0031   ~PrintMTDSens() override;
0032
0033 private:
0034   void update(const BeginOfRun *run) override;
0035   int dumpTouch(G4VPhysicalVolume *pv, unsigned int leafDepth, int ns);
0036   G4VPhysicalVolume *getTopPV();
0037   bool getBaseNumber();
0038
0039 private:
0040   std::vector<std::string> names_;
0041   std::vector<size_t> nsens_;
0042   G4NavigationHistory fHistory;
0043
0044   MTDBaseNumber thisN_;
0045   BTLNumberingScheme btlNS_;
0046   ETLNumberingScheme etlNS_;
0047 };
0048
0049 PrintMTDSens::PrintMTDSens(const edm::ParameterSet &p) : thisN_(), btlNS_(), etlNS_() {
0050   names_ = p.getUntrackedParameter<std::vector<std::string>>("Name");
0051   nsens_.resize(names_.size());
0052   for (size_t index = 0; index < nsens_.size(); index++) {
0053     nsens_[index] = 0;
0054   }
0055   G4cout << "PrintMTDSens:: Print position of MTD Sensitive Touchables: " << G4endl;
0056   for (const auto &thisName : names_) {
0057     G4cout << " for name " << thisName << "\n";
0058   }
0059   G4cout << " Total of " << names_.size() << " sensitive volume types" << G4endl;
0060 }
0061
0062 PrintMTDSens::~PrintMTDSens() {}
0063
0064 void PrintMTDSens::update(const BeginOfRun *run) {
0065   G4VPhysicalVolume *theTopPV = getTopPV();
0066   int ntotal = dumpTouch(theTopPV, 0, 0);
0067   G4cout << "\nTotal number of sensitive detector volumes for MTD is " << ntotal << G4endl;
0068   for (size_t index = 0; index < nsens_.size(); index++) {
0069     G4cout << "Sensitive volume " << names_[index] << " # copies " << nsens_[index] << G4endl;
0070   }
0071 }
0072
0073 using cms_rounding::roundIfNear0;
0074
0075 int PrintMTDSens::dumpTouch(G4VPhysicalVolume *pv, unsigned int leafDepth, int ns) {
0076   if (leafDepth == 0)
0077     fHistory.SetFirstEntry(pv);
0078   else
0079     fHistory.NewLevel(pv, kNormal, pv->GetCopyNo());
0080
0081   int ntotal(ns);
0082   G4ThreeVector globalpoint = fHistory.GetTopTransform().Inverse().TransformPoint(G4ThreeVector(0, 0, 0));
0083   G4LogicalVolume *lv = pv->GetLogicalVolume();
0084
0085   std::string mother = "World";
0086   bool printIt(false);
0087   if (pv->GetMotherLogical()) {
0088     mother = pv->GetMotherLogical()->GetName();
0089   }
0090   size_t index(0);
0091   for (const auto &thisName : names_) {
0092     G4String g4name(thisName);
0093     if (G4StrUtil::contains(lv->GetName(), g4name)) {
0094       printIt = true;
0095       break;
0096     }
0097     index++;
0098   }
0099
0100   if (lv->GetSensitiveDetector() && printIt) {
0101     std::stringstream sunitt;
0102     ++ntotal;
0103     ++nsens_[index];
0104     bool isBarrel = getBaseNumber();
0105     if (isBarrel) {
0106       BTLDetId theId(btlNS_.getUnitID(thisN_));
0107       sunitt << theId.rawId();
0108 #ifdef EDM_ML_DEBUG
0109       G4cout << theId << G4endl;
0110 #endif
0111     } else {
0112       ETLDetId theId(etlNS_.getUnitID(thisN_));
0113       sunitt << theId.rawId();
0114 #ifdef EDM_ML_DEBUG
0115       G4cout << theId << G4endl;
0116 #endif
0117     }
0118
0119     auto fround = [&](double in) {
0120       std::stringstream ss;
0121       ss << std::fixed << std::setw(14) << roundIfNear0(in);
0122       return ss.str();
0123     };
0124
0125     G4Box *thisSens = static_cast<G4Box *>(lv->GetSolid());
0126     G4ThreeVector cn1Global = fHistory.GetTopTransform().Inverse().TransformPoint(
0127         G4ThreeVector(thisSens->GetXHalfLength(), thisSens->GetYHalfLength(), thisSens->GetZHalfLength()));
0128
0129     sunitt << fround(globalpoint.x()) << fround(globalpoint.y()) << fround(globalpoint.z()) << fround(cn1Global.x())
0130            << fround(cn1Global.y()) << fround(cn1Global.z());
0131     edm::LogVerbatim("MTDG4sensUnitTest") << sunitt.str();
0132   }
0133
0134   int NoDaughters = lv->GetNoDaughters();
0135   while ((NoDaughters--) > 0) {
0136     G4VPhysicalVolume *pvD = lv->GetDaughter(NoDaughters);
0137     if (!pvD->IsReplicated())
0138       ntotal = dumpTouch(pvD, leafDepth + 1, ntotal);
0139   }
0140
0141   if (leafDepth > 0)
0142     fHistory.BackLevel();
0143   return ntotal;
0144 }
0145
0146 G4VPhysicalVolume *PrintMTDSens::getTopPV() {
0147   return G4TransportationManager::GetTransportationManager()->GetNavigatorForTracking()->GetWorldVolume();
0148 }
0149
0150 bool PrintMTDSens::getBaseNumber() {
0151   bool isBTL(false);
0152   thisN_.reset();
0153   thisN_.setSize(fHistory.GetMaxDepth());
0154   int theSize = fHistory.GetDepth() + 1;
0155   if (thisN_.getCapacity() < theSize)
0156     thisN_.setSize(theSize);
0157   //Get name and copy numbers
0158   if (theSize > 1) {
0159 #ifdef EDM_ML_DEBUG
0160     G4cout << "Building MTD basenumber:" << G4endl;
0161 #endif
0162     for (int ii = theSize; ii-- > 0;) {
0163       thisN_.addLevel(fHistory.GetVolume(ii)->GetName(), fHistory.GetReplicaNo(ii));
0164 #ifdef EDM_ML_DEBUG
0165       G4cout << "PrintMTDSens::getBaseNumber(): Adding level " << theSize - 1 - ii << ": "
0166              << fHistory.GetVolume(ii)->GetName() << "[" << fHistory.GetReplicaNo(ii) << "]" << G4endl;
0167 #endif
0168       if (!isBTL) {
0169         isBTL = G4StrUtil::contains(fHistory.GetVolume(ii)->GetName(), "BarrelTimingLayer");
0170       }
0171     }
0172   }
0173   return isBTL;
0174 }
0175
0176 #include "SimG4Core/Watcher/interface/SimWatcherFactory.h"
0177 #include "FWCore/PluginManager/interface/ModuleDef.h"
0178
0179 DEFINE_SIMWATCHER(PrintMTDSens);