Merge pull request #27588 from ianna/dd4hep-sim-dddworld-with-g4prodc…

…uts-v1

DD4hep:Simulation - DDDWorld Defines Regions and Sets G4 Production Cuts

DetectorDescription/DDCMS/interface/DDCompactView.h

@@ -21,14 +21,16 @@
 //
 //
 
-namespace cms {
+#include "DetectorDescription/DDCMS/interface/DDSpecParRegistry.h"
+#include "DetectorDescription/DDCMS/interface/DDDetector.h"
 
-  class DDDetector;
+namespace cms {
 
   class DDCompactView {
   public:
     DDCompactView(const cms::DDDetector& det) : m_det(det) {}
     const cms::DDDetector* detector() const { return &m_det; }
+    DDSpecParRegistry const& specpars() const { return m_det.specpars(); }
 
   private:
     const cms::DDDetector& m_det;

DetectorDescription/DDCMS/interface/Filter.h

@@ -34,8 +34,9 @@ namespace cms {
     int contains(std::string_view, std::string_view);
     bool isRegex(std::string_view);
     bool compareEqual(std::string_view, std::string_view);
-    std::string_view realTopName(std::string_view input);
+    std::string_view realTopName(std::string_view);
     std::vector<std::string_view> split(std::string_view, const char*);
+    std::string_view noNamespace(std::string_view);
   }  // namespace dd
 }  // namespace cms
 

DetectorDescription/DDCMS/src/Filter.cc

@@ -57,5 +57,12 @@ namespace cms {
         ret.emplace_back(str.substr(start, str.length() - start));
       return ret;
     }
+
+    std::string_view noNamespace(std::string_view input) {
+      std::string_view v = input;
+      auto first = v.find_first_of(":");
+      v.remove_prefix(std::min(first + 1, v.size()));
+      return v;
+    }
   }  // namespace dd
 }  // namespace cms

SimG4Core/Geometry/interface/DDG4ProductionCuts.h

@@ -2,6 +2,8 @@
 #define SimG4Core_DDG4ProductionCuts_H
 
 #include "SimG4Core/Geometry/interface/G4LogicalVolumeToDDLogicalPartMap.h"
+#include "DetectorDescription/DDCMS/interface/DDSpecParRegistry.h"
+#include "DDG4/Geant4GeometryInfo.h"
 
 #include <string>
 #include <vector>
@@ -12,15 +14,37 @@ class G4LogicalVolume;
 
 class DDG4ProductionCuts {
 public:
-  explicit DDG4ProductionCuts(const G4LogicalVolumeToDDLogicalPartMap&, int, bool);
+  explicit DDG4ProductionCuts(const G4LogicalVolumeToDDLogicalPartMap*, int, bool);
+
+  // ---------------------------------
+  // DD4hep specific constructor...
+  explicit DDG4ProductionCuts(const cms::DDSpecParRegistry*,
+                              const dd4hep::sim::Geant4GeometryMaps::VolumeMap*,
+                              int,
+                              bool);
+
   ~DDG4ProductionCuts();
 
 private:
   void initialize();
   void setProdCuts(const DDLogicalPart, G4Region*);
 
-  const G4LogicalVolumeToDDLogicalPartMap& map_;
+  const G4LogicalVolumeToDDLogicalPartMap* map_ = nullptr;
   G4LogicalVolumeToDDLogicalPartMap::Vector vec_;
+
+  // ---------------------------------
+  // DD4hep specific initialization,
+  //   methods, and local variables...
+  void dd4hepInitialize();
+  void setProdCuts(const cms::DDSpecPar*, G4Region*);
+
+  const dd4hep::sim::Geant4GeometryMaps::VolumeMap* dd4hepMap_ = nullptr;
+  std::vector<std::pair<G4LogicalVolume*, const cms::DDSpecPar*>> dd4hepVec_;
+  const cms::DDSpecParRegistry* specPars_;
+  cms::DDSpecParRefs specs_;
+  // ... end here.
+  // ---------------------------------
+
   const std::string keywordRegion_;
   const int verbosity_;
   const bool protonCut_;

SimG4Core/Geometry/src/DDDWorld.cc

@@ -25,7 +25,7 @@ DDDWorld::DDDWorld(const DDCompactView *pDD,
                    int verb,
                    bool cuts,
                    bool pcut) {
-  LogVerbatim("SimG4CoreApplication") << "DDDWorld: initialisation of Geant4 geometry";
+  LogVerbatim("SimG4CoreApplication") << "+++ DDDWorld: initialisation of Geant4 geometry";
   if (pDD4hep) {
     // DD4Hep
     const cms::DDDetector *det = pDD4hep->detector();
@@ -37,7 +37,9 @@ DDDWorld::DDDWorld(const DDCompactView *pDD,
     Geant4GeometryInfo *geometry = g4Geo.create(world).detach();
     lvMap = geometry->g4Volumes;
     m_world = geometry->world();
-
+    if (cuts) {
+      DDG4ProductionCuts pcuts(&det->specpars(), &lvMap, verb, pcut);
+    }
   } else {
     // old DD code
     G4LogicalVolumeToDDLogicalPartMap lvMap;
@@ -47,7 +49,7 @@ DDDWorld::DDDWorld(const DDCompactView *pDD,
     LogVerbatim("SimG4CoreApplication") << "DDDWorld: worldLV: " << world->GetName();
     m_world = new G4PVPlacement(nullptr, G4ThreeVector(), world, "DDDWorld", nullptr, false, 0);
     if (cuts) {
-      DDG4ProductionCuts pcuts(lvMap, verb, pcut);
+      DDG4ProductionCuts pcuts(&lvMap, verb, pcut);
     }
   }
   LogVerbatim("SimG4CoreApplication") << "DDDWorld: initialisation of Geant4 geometry is done.";

SimG4Core/Geometry/src/DDG4ProductionCuts.cc

@@ -1,5 +1,6 @@
 #include "SimG4Core/Geometry/interface/DDG4ProductionCuts.h"
 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
+#include "DetectorDescription/DDCMS/interface/Filter.h"
 
 #include "FWCore/MessageLogger/interface/MessageLogger.h"
 #include "FWCore/Utilities/interface/Exception.h"
@@ -11,39 +12,58 @@
 
 #include <algorithm>
 
-DDG4ProductionCuts::DDG4ProductionCuts(const G4LogicalVolumeToDDLogicalPartMap& map, int verb, bool pcut)
-    : map_(map), keywordRegion_("CMSCutsRegion"), verbosity_(verb), protonCut_(pcut) {
-  initialize();
-}
-
-DDG4ProductionCuts::~DDG4ProductionCuts() {}
-
-/** helper function to compare parts through their name instead of comparing them
-    by their pointers. 
-    It's guaranteed to produce the same order in subsequent application runs,
-    while pointers usually can't guarantee this
-*/
-bool dd_is_greater(const std::pair<G4LogicalVolume*, DDLogicalPart>& p1,
-                   const std::pair<G4LogicalVolume*, DDLogicalPart>& p2) {
-  bool result = false;
-  if (p1.second.name().ns() > p2.second.name().ns()) {
-    result = true;
-  }
-  if (p1.second.name().ns() == p2.second.name().ns()) {
-    if (p1.second.name().name() > p2.second.name().name()) {
+namespace {
+  /** helper function to compare parts through their name instead of comparing them
+      by their pointers. 
+      It's guaranteed to produce the same order in subsequent application runs,
+      while pointers usually can't guarantee this
+  */
+  bool dd_is_greater(const std::pair<G4LogicalVolume*, DDLogicalPart>& p1,
+                     const std::pair<G4LogicalVolume*, DDLogicalPart>& p2) {
+    bool result = false;
+    if (p1.second.name().ns() > p2.second.name().ns()) {
       result = true;
     }
-    if (p1.second.name().name() == p2.second.name().name()) {
-      if (p1.first->GetName() > p2.first->GetName()) {
+    if (p1.second.name().ns() == p2.second.name().ns()) {
+      if (p1.second.name().name() > p2.second.name().name()) {
         result = true;
       }
+      if (p1.second.name().name() == p2.second.name().name()) {
+        if (p1.first->GetName() > p2.first->GetName()) {
+          result = true;
+        }
+      }
     }
+    return result;
   }
-  return result;
+
+  bool sortByName(const std::pair<G4LogicalVolume*, const cms::DDSpecPar*>& p1,
+                  const std::pair<G4LogicalVolume*, const cms::DDSpecPar*>& p2) {
+    bool result = false;
+    if (p1.first->GetName() > p2.first->GetName()) {
+      result = true;
+    }
+    return result;
+  }
+}  // namespace
+
+DDG4ProductionCuts::DDG4ProductionCuts(const G4LogicalVolumeToDDLogicalPartMap* map, int verb, bool pcut)
+    : map_(map), keywordRegion_("CMSCutsRegion"), verbosity_(verb), protonCut_(pcut) {
+  initialize();
 }
 
+DDG4ProductionCuts::DDG4ProductionCuts(const cms::DDSpecParRegistry* specPars,
+                                       const dd4hep::sim::Geant4GeometryMaps::VolumeMap* map,
+                                       int verb,
+                                       bool pcut)
+    : dd4hepMap_(map), specPars_(specPars), keywordRegion_("CMSCutsRegion"), verbosity_(verb), protonCut_(pcut) {
+  dd4hepInitialize();
+}
+
+DDG4ProductionCuts::~DDG4ProductionCuts() {}
+
 void DDG4ProductionCuts::initialize() {
-  vec_ = map_.all(keywordRegion_);
+  vec_ = map_->all(keywordRegion_);
   // sort all root volumes - to get the same sequence at every run of the application.
   // (otherwise, the sequence will depend on the pointer (memory address) of the
   // involved objects, because 'new' does no guarantee that you allways get a
@@ -62,11 +82,11 @@ void DDG4ProductionCuts::initialize() {
   G4Region* region = nullptr;
   G4RegionStore* store = G4RegionStore::GetInstance();
   for (auto const& vv : vec_) {
-    unsigned int num = map_.toString(keywordRegion_, vv.second, regionName);
+    unsigned int num = map_->toString(keywordRegion_, vv.second, regionName);
     edm::LogVerbatim("Geometry") << "  num  " << num << " regionName: " << regionName << " " << store;
 
     if (num != 1) {
-      throw cms::Exception("SimG4CorePhysics", " DDG4ProductionCuts::initialize: Problem with Region tags.");
+      throw cms::Exception("SimG4CoreGeometry", " DDG4ProductionCuts::initialize: Problem with Region tags.");
     }
     if (regionName != curName) {
       edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : regionName " << regionName << " " << store;
@@ -87,6 +107,45 @@ void DDG4ProductionCuts::initialize() {
   }
 }
 
+void DDG4ProductionCuts::dd4hepInitialize() {
+  specPars_->filter(specs_, keywordRegion_);
+
+  for (auto const& it : *dd4hepMap_) {
+    for (auto const& fit : specs_) {
+      for (auto const& pit : fit->paths) {
+        if (cms::dd::compareEqual(cms::dd::noNamespace(it.first.name()), cms::dd::realTopName(pit))) {
+          dd4hepVec_.emplace_back(std::make_pair<G4LogicalVolume*, const cms::DDSpecPar*>(&*it.second, &*fit));
+        }
+      }
+    }
+  }
+  // sort all root volumes - to get the same sequence at every run of the application.
+  sort(begin(dd4hepVec_), end(dd4hepVec_), &sortByName);
+
+  // Now generate all the regions
+  for (auto const& it : dd4hepVec_) {
+    auto regName = it.second->strValue(keywordRegion_.data());
+    G4Region* region = G4RegionStore::GetInstance()->FindOrCreateRegion({regName.data(), regName.size()});
+    region->AddRootLogicalVolume(it.first);
+    edm::LogVerbatim("Geometry") << it.first->GetName() << ": " << it.second->strValue(keywordRegion_.data());
+    edm::LogVerbatim("Geometry") << " MakeRegions: added " << it.first->GetName() << " to region " << region->GetName();
+    edm::LogVerbatim("Geometry").log([&](auto& log) {
+      for (auto const& sit : it.second->spars) {
+        log << sit.first << " =  " << sit.second[0] << "\n";
+      }
+    });
+    setProdCuts(it.second, region);
+  }
+
+  if (verbosity_ > 0) {
+    LogDebug("Geometry") << " DDG4ProductionCuts (New) : starting\n"
+                         << " DDG4ProductionCuts : Got " << dd4hepVec_.size() << " region roots.\n"
+                         << " DDG4ProductionCuts : List of all roots:";
+    for (size_t jj = 0; jj < dd4hepVec_.size(); ++jj)
+      LogDebug("Geometry") << "   DDG4ProductionCuts : root=" << dd4hepVec_[jj];
+  }
+}
+
 void DDG4ProductionCuts::setProdCuts(const DDLogicalPart lpart, G4Region* region) {
   //
   // search for production cuts
@@ -96,25 +155,25 @@ void DDG4ProductionCuts::setProdCuts(const DDLogicalPart lpart, G4Region* region
   double electroncut = 0.0;
   double positroncut = 0.0;
   double protoncut = 0.0;
-  int temp = map_.toDouble("ProdCutsForGamma", lpart, gammacut);
+  int temp = map_->toDouble("ProdCutsForGamma", lpart, gammacut);
   if (temp != 1) {
     throw cms::Exception(
         "SimG4CorePhysics",
         " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForGamma.");
   }
-  temp = map_.toDouble("ProdCutsForElectrons", lpart, electroncut);
+  temp = map_->toDouble("ProdCutsForElectrons", lpart, electroncut);
   if (temp != 1) {
     throw cms::Exception(
         "SimG4CorePhysics",
         " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForElectrons.");
   }
-  temp = map_.toDouble("ProdCutsForPositrons", lpart, positroncut);
+  temp = map_->toDouble("ProdCutsForPositrons", lpart, positroncut);
   if (temp != 1) {
     throw cms::Exception(
         "SimG4CorePhysics",
         " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForPositrons.");
   }
-  temp = map_.toDouble("ProdCutsForProtons", lpart, protoncut);
+  temp = map_->toDouble("ProdCutsForProtons", lpart, protoncut);
   if (temp == 0) {
     // There is no ProdCutsForProtons set in XML,
     // check if it's a legacy geometry scenario without it
@@ -147,3 +206,87 @@ void DDG4ProductionCuts::setProdCuts(const DDLogicalPart lpart, G4Region* region
                                  << "\n    Gamma    : " << gammacut << "\n    Proton   : " << protoncut;
   }
 }
+
+void DDG4ProductionCuts::setProdCuts(const cms::DDSpecPar* spec, G4Region* region) {
+  //
+  // Create and fill production cuts
+  //
+  G4ProductionCuts* prodCuts = region->GetProductionCuts();
+  if (!prodCuts) {
+    // FIXME: Here we use a dd4hep string to double evaluator
+    //        Beware of the units!!!
+
+    //
+    // search for production cuts
+    // you must have four of them: e+ e- gamma proton
+    //
+    double gammacut = 0.0;
+    double electroncut = 0.0;
+    double positroncut = 0.0;
+    double protoncut = 0.0;
+
+    auto gammacutStr = spec->strValue("ProdCutsForGamma");
+    if (gammacutStr.empty()) {
+      throw cms::Exception(
+          "SimG4CorePhysics",
+          " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForGamma.");
+    }
+    gammacut = dd4hep::_toDouble({gammacutStr.data(), gammacutStr.size()});
+
+    auto electroncutStr = spec->strValue("ProdCutsForElectrons");
+    if (electroncutStr.empty()) {
+      throw cms::Exception(
+          "SimG4CorePhysics",
+          " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForElectrons.");
+    }
+    electroncut = dd4hep::_toDouble({electroncutStr.data(), electroncutStr.size()});
+
+    auto positroncutStr = spec->strValue("ProdCutsForPositrons");
+    if (positroncutStr.empty()) {
+      throw cms::Exception(
+          "SimG4CorePhysics",
+          " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForPositrons.");
+    }
+    positroncut = dd4hep::_toDouble({positroncutStr.data(), positroncutStr.size()});
+
+    if (!spec->hasValue("ProdCutsForProtons")) {
+      // There is no ProdCutsForProtons set in XML,
+      // check if it's a legacy geometry scenario without it
+      if (protonCut_) {
+        protoncut = electroncut;
+      } else {
+        protoncut = 0.;
+      }
+    } else {
+      auto protoncutStr = spec->strValue("ProdCutsForProtons");
+      if (protoncutStr.empty()) {
+        throw cms::Exception(
+            "SimG4CorePhysics",
+            " DDG4ProductionCuts::setProdCuts: Problem with Region tags - more than one ProdCutsForProtons.");
+      }
+      protoncut = dd4hep::_toDouble({protoncutStr.data(), protoncutStr.size()});
+    }
+
+    prodCuts = new G4ProductionCuts();
+    region->SetProductionCuts(prodCuts);
+
+    prodCuts->SetProductionCut(gammacut, idxG4GammaCut);
+    prodCuts->SetProductionCut(electroncut, idxG4ElectronCut);
+    prodCuts->SetProductionCut(positroncut, idxG4PositronCut);
+    prodCuts->SetProductionCut(protoncut, idxG4ProtonCut);
+    if (verbosity_ > 0) {
+      edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : Setting cuts for " << region->GetName()
+                                   << "\n    Electrons: " << electroncut << "\n    Positrons: " << positroncut
+                                   << "\n    Gamma    : " << gammacut << "\n    Proton   : " << protoncut;
+    }
+  } else {
+    if (verbosity_ > 0) {
+      edm::LogVerbatim("Geometry")
+          << "DDG4ProductionCuts : Cuts are already set for " << region->GetName()
+          << "\n    Electrons: " << region->GetProductionCuts()->GetProductionCut(idxG4ElectronCut)
+          << "\n    Positrons: " << region->GetProductionCuts()->GetProductionCut(idxG4PositronCut)
+          << "\n    Gamma    : " << region->GetProductionCuts()->GetProductionCut(idxG4GammaCut)
+          << "\n    Proton   : " << region->GetProductionCuts()->GetProductionCut(idxG4ProtonCut);
+    }
+  }
+}