test/plugins/DD4hepTestDDDWorld.cc

0001 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0002 #include "FWCore/Framework/interface/MakerMacros.h"
0003 #include "FWCore/Framework/interface/ESTransientHandle.h"
0004 #include "FWCore/Framework/interface/EventSetup.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "DataFormats/Math/interface/GeantUnits.h"
0007 #include "DetectorDescription/DDCMS/interface/DDDetector.h"
0008 #include "DetectorDescription/DDCMS/interface/DDVectorRegistry.h"
0009 #include "Geometry/Records/interface/DDVectorRegistryRcd.h"
0010 #include "Geometry/Records/interface/IdealGeometryRecord.h"
0011 #include "Geometry/Records/interface/DDSpecParRegistryRcd.h"
0012 #include "SimG4Core/DD4hepGeometry/interface/DD4hep_DDDWorld.h"
0013 #include <DDG4/Geant4Converter.h>
0014 #include <DD4hep/Detector.h>
0015 #include <DD4hep/Handle.h>
0016 #include <DD4hep/Filter.h>
0017 #include <DD4hep/SpecParRegistry.h>
0018 #include "G4LogicalVolume.hh"
0019 #include "G4MTRunManagerKernel.hh"
0020 #include "G4ProductionCuts.hh"
0021 #include "G4RegionStore.hh"
0022
0023 #include <iostream>
0024 #include <string>
0025 #include <string_view>
0026
0027 using namespace std;
0028 using namespace cms;
0029 using namespace edm;
0030 using namespace dd4hep;
0031
0032 namespace {
0033   bool sortByName(const std::pair<G4LogicalVolume*, const dd4hep::SpecPar*>& p1,
0034                   const std::pair<G4LogicalVolume*, const dd4hep::SpecPar*>& p2) {
0035     bool result = false;
0036     if (p1.first->GetName() > p2.first->GetName()) {
0037       result = true;
0038     }
0039     return result;
0040   }
0041 }  // namespace
0042
0043 class DD4hepTestDDDWorld : public one::EDAnalyzer<> {
0044 public:
0045   explicit DD4hepTestDDDWorld(const ParameterSet&);
0046
0047   void beginJob() override {}
0048   void analyze(Event const& iEvent, EventSetup const&) override;
0049   void endJob() override {}
0050
0051 private:
0052   void update();
0053   void initialize(const dd4hep::sim::Geant4GeometryMaps::VolumeMap&);
0054   G4ProductionCuts* getProductionCuts(G4Region* region);
0055
0056   const ESInputTag tag_;
0057   const ESGetToken<DDVectorRegistry, DDVectorRegistryRcd> vecRegToken_;
0058   const ESGetToken<DDDetector, IdealGeometryRecord> detectorToken_;
0059   const ESGetToken<dd4hep::SpecParRegistry, DDSpecParRegistryRcd> registryToken_;
0060   const dd4hep::SpecParRegistry* specPars_;
0061   G4MTRunManagerKernel* kernel_;
0062   dd4hep::SpecParRefs specs_;
0063   vector<pair<G4LogicalVolume*, const dd4hep::SpecPar*>> vec_;
0064   const string keywordRegion_;
0065   unique_ptr<DDDWorld> world_;
0066   int verbosity_;
0067 };
0068
0069 DD4hepTestDDDWorld::DD4hepTestDDDWorld(const ParameterSet& iConfig)
0070     : tag_(iConfig.getParameter<ESInputTag>("DDDetector")),
0071       vecRegToken_(esConsumes(tag_)),
0072       detectorToken_(esConsumes(tag_)),
0073       registryToken_(esConsumes(tag_)),
0074       keywordRegion_("CMSCutsRegion") {
0075   verbosity_ = iConfig.getUntrackedParameter<int>("Verbosity", 1);
0076   kernel_ = new G4MTRunManagerKernel();
0077 }
0078
0079 void DD4hepTestDDDWorld::analyze(const Event&, const EventSetup& iEventSetup) {
0080   LogVerbatim("Geometry") << "\nDD4hepTestDDDWorld::analyze: " << tag_;
0081
0082   ESTransientHandle<DDVectorRegistry> reg = iEventSetup.getTransientHandle(vecRegToken_);
0083
0084   ESTransientHandle<DDDetector> ddd = iEventSetup.getTransientHandle(detectorToken_);
0085
0086   ESTransientHandle<dd4hep::SpecParRegistry> registry = iEventSetup.getTransientHandle(registryToken_);
0087   specPars_ = registry.product();
0088
0089   const dd4hep::Detector& detector = *ddd->description();
0090   dd4hep::sim::Geant4Converter g4Geo = dd4hep::sim::Geant4Converter(detector);
0091   g4Geo.debugMaterials = true;
0092   g4Geo.debugElements = true;
0093   g4Geo.debugShapes = true;
0094   g4Geo.debugVolumes = true;
0095   g4Geo.debugPlacements = true;
0096   g4Geo.create(detector.world());
0097
0098   dd4hep::sim::Geant4GeometryMaps::VolumeMap lvMap;
0099   world_.reset(new DDDWorld(ddd.product(), lvMap));
0100   initialize(lvMap);
0101   update();
0102   LogVerbatim("Geometry") << "Done.";
0103 }
0104
0105 void DD4hepTestDDDWorld::initialize(const dd4hep::sim::Geant4GeometryMaps::VolumeMap& vmap) {
0106   specPars_->filter(specs_, keywordRegion_);
0107
0108   LogVerbatim("Geometry").log([&](auto& log) {
0109     for (auto const& it : vmap) {
0110       for (auto const& fit : specs_) {
0111         for (auto const& sit : fit.second->numpars) {
0112           log << sit.first << " =  " << sit.second[0] << "\n";
0113         }
0114         for (auto const& pit : fit.second->paths) {
0115           const std::string_view selection = dd4hep::dd::realTopName(pit);
0116           const std::string_view name = dd4hep::dd::noNamespace(it.first.name());
0117           log << selection << "\n";
0118           log << "   compare equal to " << name << " ... ";
0119           if (!(dd4hep::dd::isRegex(selection))
0120                   ? dd4hep::dd::compareEqual(name, selection)
0121                   : regex_match(name.begin(), name.end(), regex(std::string(selection)))) {
0122             vec_.emplace_back(std::make_pair<G4LogicalVolume*, const dd4hep::SpecPar*>(&*it.second, &*fit.second));
0123             log << "   are equal!\n";
0124           } else
0125             log << "   nope.\n";
0126         }
0127       }
0128     }
0129   });
0130
0131   // sort all root volumes - to get the same sequence at every run of the application.
0132   sort(begin(vec_), end(vec_), &sortByName);
0133
0134   // Now generate all the regions
0135   for (auto const& it : vec_) {
0136     auto regName = it.second->strValue(keywordRegion_);
0137     G4Region* region = G4RegionStore::GetInstance()->FindOrCreateRegion({regName.data(), regName.size()});
0138     region->AddRootLogicalVolume(it.first);
0139     LogVerbatim("Geometry") << it.first->GetName() << ": " << it.second->strValue(keywordRegion_);
0140     LogVerbatim("Geometry") << " MakeRegions: added " << it.first->GetName() << " to region " << region->GetName();
0141   }
0142 }
0143
0144 void DD4hepTestDDDWorld::update() {
0145   LogVerbatim("Geometry").log([&](auto& log) {
0146     log << "DD4hepTestDDDWorld::update()\n";
0147     for (const auto& t : vec_) {
0148       log << t.first->GetName() << ":\n";
0149       for (const auto& kl : t.second->numpars) {
0150         log << kl.first << " = ";
0151         for (const auto& kil : kl.second) {
0152           log << kil << ", ";
0153         }
0154       }
0155       log << "\n";
0156     }
0157     log << "DD4hepTestDDDWorld::update() done!\n";
0158   });
0159   // Loop over all DDLP and provide the cuts for each region
0160   for (auto const& it : vec_) {
0161     auto regName = it.second->strValue(keywordRegion_);
0162     G4Region* region = G4RegionStore::GetInstance()->FindOrCreateRegion({regName.data(), regName.size()});
0163
0164     //
0165     // search for production cuts
0166     // you must have four of them: e+ e- gamma proton
0167     //
0168
0169     // Geant4 expects mm units. DD4hep may return different units, so convert to mm.
0170     double gammacut = it.second->dblValue("ProdCutsForGamma") / dd4hep::mm;
0171     double electroncut = it.second->dblValue("ProdCutsForElectrons") / dd4hep::mm;
0172     double positroncut = it.second->dblValue("ProdCutsForPositrons") / dd4hep::mm;
0173     double protoncut = it.second->dblValue("ProdCutsForProtons") / dd4hep::mm;
0174     if (protoncut == 0) {
0175       protoncut = electroncut;
0176     }
0177
0178     //
0179     // For the moment I assume all of the four are set
0180     //
0181     G4ProductionCuts* prodCuts = getProductionCuts(region);
0182     prodCuts->SetProductionCut(gammacut, idxG4GammaCut);
0183     prodCuts->SetProductionCut(electroncut, idxG4ElectronCut);
0184     prodCuts->SetProductionCut(positroncut, idxG4PositronCut);
0185     prodCuts->SetProductionCut(protoncut, idxG4ProtonCut);
0186     if (verbosity_ > 0) {
0187       LogVerbatim("Geometry") << "DDG4ProductionCuts : Setting cuts for " << regName
0188                               << "\n    Electrons: " << electroncut << " mm\n    Positrons: " << positroncut
0189                               << " mm\n    Gamma    : " << gammacut << " mm\n    Protons  : " << protoncut << " mm\n";
0190     }
0191   }
0192 }
0193
0194 G4ProductionCuts* DD4hepTestDDDWorld::getProductionCuts(G4Region* region) {
0195   G4ProductionCuts* prodCuts = region->GetProductionCuts();
0196   if (!prodCuts) {
0197     prodCuts = new G4ProductionCuts();
0198     region->SetProductionCuts(prodCuts);
0199   }
0200   return prodCuts;
0201 }
0202
0203 DEFINE_FWK_MODULE(DD4hepTestDDDWorld);