Geometry/src/DDG4ProductionCuts.cc

0001 #include "SimG4Core/Geometry/interface/DDG4ProductionCuts.h"
0002 #include "DetectorDescription/Core/interface/DDLogicalPart.h"
0003
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/Utilities/interface/Exception.h"
0006
0007 #include <DD4hep/Filter.h>
0008 #include <DD4hep/SpecParRegistry.h>
0009
0010 #include "G4ProductionCuts.hh"
0011 #include "G4RegionStore.hh"
0012 #include "G4Region.hh"
0013 #include "G4LogicalVolume.hh"
0014 #include "G4LogicalVolumeStore.hh"
0015
0016 #include <algorithm>
0017
0018 namespace {
0019   /** helper function to compare parts through their name instead of comparing them
0020       by their pointers.
0021       It's guaranteed to produce the same order in subsequent application runs,
0022       while pointers usually can't guarantee this
0023   */
0024   bool dd_is_greater(const std::pair<G4LogicalVolume*, DDLogicalPart>& p1,
0025                      const std::pair<G4LogicalVolume*, DDLogicalPart>& p2) {
0026     bool result = false;
0027     if (p1.second.name().ns() > p2.second.name().ns()) {
0028       result = true;
0029     }
0030     if (p1.second.name().ns() == p2.second.name().ns()) {
0031       if (p1.second.name().name() > p2.second.name().name()) {
0032         result = true;
0033       }
0034       if (p1.second.name().name() == p2.second.name().name()) {
0035         if (p1.first->GetName() > p2.first->GetName()) {
0036           result = true;
0037         }
0038       }
0039     }
0040     return result;
0041   }
0042
0043   bool sortByName(const std::pair<G4LogicalVolume*, const dd4hep::SpecPar*>& p1,
0044                   const std::pair<G4LogicalVolume*, const dd4hep::SpecPar*>& p2) {
0045     bool result = false;
0046     if (p1.first->GetName() > p2.first->GetName()) {
0047       result = true;
0048     }
0049     return result;
0050   }
0051 }  // namespace
0052
0053 DDG4ProductionCuts::DDG4ProductionCuts(const G4LogicalVolumeToDDLogicalPartMap* map, int verb, bool pcut)
0054     : map_(map), keywordRegion_("CMSCutsRegion"), verbosity_(verb), protonCut_(pcut) {
0055   initialize();
0056 }
0057
0058 DDG4ProductionCuts::DDG4ProductionCuts(const dd4hep::SpecParRegistry* specPars,
0059                                        const dd4hep::sim::Geant4GeometryMaps::VolumeMap* map,
0060                                        int verb,
0061                                        bool pcut)
0062     : dd4hepMap_(map), specPars_(specPars), keywordRegion_("CMSCutsRegion"), verbosity_(verb), protonCut_(pcut) {
0063   dd4hepInitialize();
0064 }
0065
0066 void DDG4ProductionCuts::initialize() {
0067   vec_ = map_->all(keywordRegion_);
0068   // sort all root volumes - to get the same sequence at every run of the application.
0069   // (otherwise, the sequence will depend on the pointer (memory address) of the
0070   // involved objects, because 'new' does no guarantee that you allways get a
0071   // higher (or lower) address when allocating an object of the same type ...
0072   sort(vec_.begin(), vec_.end(), &dd_is_greater);
0073   if (verbosity_ > 0) {
0074     edm::LogVerbatim("Geometry") << " DDG4ProductionCuts : got " << vec_.size() << " region roots.\n"
0075                                  << " DDG4ProductionCuts : List of all roots:";
0076     for (auto const& vv : vec_)
0077       edm::LogVerbatim("Geometry") << "    " << vv.first->GetName() << " : " << vv.second.name();
0078   }
0079
0080   // Now generate all the regions
0081   std::string curName = "";
0082   std::string regionName = "";
0083   G4Region* region = nullptr;
0084   G4RegionStore* store = G4RegionStore::GetInstance();
0085   for (auto const& vv : vec_) {
0086     unsigned int num = map_->toString(keywordRegion_, vv.second, regionName);
0087     edm::LogVerbatim("Geometry") << "  num  " << num << " regionName: " << regionName << ", the store of size "
0088                                  << store->size();
0089     if (num != 1) {
0090       throw cms::Exception("SimG4CoreGeometry", " DDG4ProductionCuts::initialize: Problem with Region tags.");
0091       return;
0092     }
0093     if (regionName != curName) {
0094       edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : regionName " << regionName << ", the store of size "
0095                                    << store->size();
0096       region = store->FindOrCreateRegion(regionName);
0097       edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : region " << region->GetName();
0098       if (nullptr == region) {
0099         throw cms::Exception("SimG4CoreGeometry", " DDG4ProductionCuts::initialize: Problem with Region tags.");
0100         return;
0101       }
0102       curName = regionName;
0103       edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : new G4Region " << vv.first->GetName();
0104       setProdCuts(vv.second, region);
0105     }
0106     if (nullptr != region) {
0107       region->AddRootLogicalVolume(vv.first);
0108       if (verbosity_ > 0)
0109         edm::LogVerbatim("Geometry") << "  added " << vv.first->GetName() << " to region " << region->GetName();
0110     }
0111   }
0112 }
0113
0114 void DDG4ProductionCuts::dd4hepInitialize() {
0115   dd4hep::SpecParRefs specs;
0116   specPars_->filter(specs, keywordRegion_);
0117
0118   // LOOP ON ALL LOGICAL VOLUMES
0119   for (auto const& it : *dd4hepMap_) {
0120     bool foundMatch = false;  // Same behavior as in DDD: when matching SpecPar is found, stop search!
0121     // SEARCH ON ALL SPECPARS
0122     for (auto const& fit : specs) {
0123       for (auto const& pit : fit.second->paths) {
0124         const std::string_view selection = dd4hep::dd::noNamespace(dd4hep::dd::realTopName(pit));
0125         const std::string_view name = dd4hep::dd::noNamespace(it.first.name());
0126         if (!(dd4hep::dd::isRegex(selection))
0127                 ? dd4hep::dd::compareEqual(name, selection)
0128                 : std::regex_match(name.begin(), name.end(), std::regex(selection.begin(), selection.end()))) {
0129           dd4hepVec_.emplace_back(std::make_pair<G4LogicalVolume*, const dd4hep::SpecPar*>(&*it.second, &*fit.second));
0130           foundMatch = true;
0131           break;
0132         }
0133       }
0134       if (foundMatch)
0135         break;
0136     }  // Search on all SpecPars
0137   }    // Loop on all logical volumes
0138
0139   // sort all root volumes - to get the same sequence at every run of the application.
0140   sort(begin(dd4hepVec_), end(dd4hepVec_), &sortByName);
0141
0142   // Now generate all the regions
0143   for (auto const& it : dd4hepVec_) {
0144     auto regName = it.second->strValue(keywordRegion_);
0145     G4Region* region = G4RegionStore::GetInstance()->FindOrCreateRegion({regName.data(), regName.size()});
0146
0147     region->AddRootLogicalVolume(it.first);
0148     edm::LogVerbatim("Geometry") << it.first->GetName() << ": " << regName;
0149     edm::LogVerbatim("Geometry") << " MakeRegions: added " << it.first->GetName() << " to region " << region->GetName();
0150
0151     // Also treat reflected volumes
0152     const G4String& nonReflectedG4Name = it.first->GetName();
0153     const G4String& reflectedG4Name = nonReflectedG4Name + "_refl";
0154     const G4LogicalVolumeStore* const allG4LogicalVolumes = G4LogicalVolumeStore::GetInstance();
0155     const auto reflectedG4LogicalVolumeIt = std::find_if(
0156         allG4LogicalVolumes->begin(), allG4LogicalVolumes->end(), [&](const G4LogicalVolume* const aG4LogicalVolume) {
0157           return (aG4LogicalVolume->GetName() == reflectedG4Name);
0158         });
0159     // If G4 Logical volume has a reflected volume, add it to the region as well.
0160     if (reflectedG4LogicalVolumeIt != allG4LogicalVolumes->end()) {
0161       region->AddRootLogicalVolume(*reflectedG4LogicalVolumeIt);
0162       edm::LogVerbatim("Geometry") << " MakeRegions: added " << (*reflectedG4LogicalVolumeIt)->GetName()
0163                                    << " to region " << region->GetName();
0164     }
0165
0166     edm::LogVerbatim("Geometry").log([&](auto& log) {
0167       for (auto const& sit : it.second->spars) {
0168         log << sit.first << " =  " << sit.second[0] << "\n";
0169       }
0170     });
0171     setProdCuts(it.second, region);
0172   }
0173
0174   if (verbosity_ > 0) {
0175     edm::LogVerbatim("SimG4CoreGeometry") << " DDG4ProductionCuts (New) : starting\n"
0176                                           << " DDG4ProductionCuts : Got " << dd4hepVec_.size() << " region roots.\n"
0177                                           << " DDG4ProductionCuts : List of all roots:";
0178     for (size_t jj = 0; jj < dd4hepVec_.size(); ++jj)
0179       edm::LogVerbatim("SimG4CoreGeometry") << "   DDG4ProductionCuts : root=" << dd4hepVec_[jj].first->GetName()
0180                                             << " , " << dd4hepVec_[jj].second->paths.at(0);
0181   }
0182 }
0183
0184 void DDG4ProductionCuts::setProdCuts(const DDLogicalPart lpart, G4Region* region) {
0185   //
0186   // search for production cuts
0187   // you must have four of them: e+ e- gamma proton
0188   //
0189   double gammacut = 0.0;
0190   double electroncut = 0.0;
0191   double positroncut = 0.0;
0192   double protoncut = 0.0;
0193   int temp = map_->toDouble("ProdCutsForGamma", lpart, gammacut);
0194   if (temp != 1) {
0195     throw cms::Exception(
0196         "SimG4CorePhysics",
0197         " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForGamma.");
0198   }
0199   temp = map_->toDouble("ProdCutsForElectrons", lpart, electroncut);
0200   if (temp != 1) {
0201     throw cms::Exception(
0202         "SimG4CorePhysics",
0203         " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForElectrons.");
0204   }
0205   temp = map_->toDouble("ProdCutsForPositrons", lpart, positroncut);
0206   if (temp != 1) {
0207     throw cms::Exception(
0208         "SimG4CorePhysics",
0209         " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForPositrons.");
0210   }
0211   temp = map_->toDouble("ProdCutsForProtons", lpart, protoncut);
0212   if (temp == 0) {
0213     // There is no ProdCutsForProtons set in XML,
0214     // check if it's a legacy geometry scenario without it
0215     if (protonCut_) {
0216       protoncut = electroncut;
0217     } else {
0218       protoncut = 0.;
0219     }
0220   } else if (temp != 1) {
0221     throw cms::Exception(
0222         "SimG4CorePhysics",
0223         " DDG4ProductionCuts::setProdCuts: Problem with Region tags - more than one ProdCutsForProtons.");
0224   }
0225
0226   //
0227   // Create and fill production cuts
0228   //
0229   G4ProductionCuts* prodCuts = region->GetProductionCuts();
0230   if (nullptr == prodCuts) {
0231     prodCuts = new G4ProductionCuts();
0232     region->SetProductionCuts(prodCuts);
0233   }
0234   prodCuts->SetProductionCut(gammacut, idxG4GammaCut);
0235   prodCuts->SetProductionCut(electroncut, idxG4ElectronCut);
0236   prodCuts->SetProductionCut(positroncut, idxG4PositronCut);
0237   prodCuts->SetProductionCut(protoncut, idxG4ProtonCut);
0238   if (verbosity_ > 0) {
0239     edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : Setting cuts for " << region->GetName()
0240                                  << "\n    Electrons: " << electroncut << "\n    Positrons: " << positroncut
0241                                  << "\n    Gamma    : " << gammacut << "\n    Proton   : " << protoncut;
0242   }
0243 }
0244
0245 void DDG4ProductionCuts::setProdCuts(const dd4hep::SpecPar* spec, G4Region* region) {
0246   //
0247   // Create and fill production cuts
0248   //
0249   G4ProductionCuts* prodCuts = region->GetProductionCuts();
0250   if (!prodCuts) {
0251     //
0252     // search for production cuts
0253     // you must have four of them: e+ e- gamma proton
0254     //
0255     double gammacut = spec->dblValue("ProdCutsForGamma") / dd4hep::mm;  // Convert from DD4hep units to mm
0256     double electroncut = spec->dblValue("ProdCutsForElectrons") / dd4hep::mm;
0257     double positroncut = spec->dblValue("ProdCutsForPositrons") / dd4hep::mm;
0258     double protoncut = spec->dblValue("ProdCutsForProtons") / dd4hep::mm;
0259     if (protoncut == 0) {
0260       protoncut = electroncut;
0261     }
0262
0263     prodCuts = new G4ProductionCuts();
0264     region->SetProductionCuts(prodCuts);
0265
0266     prodCuts->SetProductionCut(gammacut, idxG4GammaCut);
0267     prodCuts->SetProductionCut(electroncut, idxG4ElectronCut);
0268     prodCuts->SetProductionCut(positroncut, idxG4PositronCut);
0269     prodCuts->SetProductionCut(protoncut, idxG4ProtonCut);
0270     if (verbosity_ > 0) {
0271       edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : Setting cuts for " << region->GetName()
0272                                    << "\n    Electrons: " << electroncut << "\n    Positrons: " << positroncut
0273                                    << "\n    Gamma    : " << gammacut << "\n    Proton   : " << protoncut;
0274     }
0275   } else {
0276     if (verbosity_ > 0) {
0277       edm::LogVerbatim("Geometry")
0278           << "DDG4ProductionCuts : Cuts are already set for " << region->GetName()
0279           << "\n    Electrons: " << region->GetProductionCuts()->GetProductionCut(idxG4ElectronCut)
0280           << "\n    Positrons: " << region->GetProductionCuts()->GetProductionCut(idxG4PositronCut)
0281           << "\n    Gamma    : " << region->GetProductionCuts()->GetProductionCut(idxG4GammaCut)
0282           << "\n    Proton   : " << region->GetProductionCuts()->GetProductionCut(idxG4ProtonCut);
0283     }
0284   }
0285 }