CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
DDG4ProductionCuts.cc
Go to the documentation of this file.
4 
7 
8 #include <DD4hep/Filter.h>
9 #include <DD4hep/SpecParRegistry.h>
10 
11 #include "G4ProductionCuts.hh"
12 #include "G4RegionStore.hh"
13 #include "G4Region.hh"
14 #include "G4LogicalVolume.hh"
15 
16 #include <algorithm>
17 
19 
20 #ifdef HAVE_GEANT4_UNITS
21 #define MM_2_CM 1.0
22 #else
23 #define MM_2_CM 0.1
24 #endif
25 
26 namespace {
32  bool dd_is_greater(const std::pair<G4LogicalVolume*, DDLogicalPart>& p1,
33  const std::pair<G4LogicalVolume*, DDLogicalPart>& p2) {
34  bool result = false;
35  if (p1.second.name().ns() > p2.second.name().ns()) {
36  result = true;
37  }
38  if (p1.second.name().ns() == p2.second.name().ns()) {
39  if (p1.second.name().name() > p2.second.name().name()) {
40  result = true;
41  }
42  if (p1.second.name().name() == p2.second.name().name()) {
43  if (p1.first->GetName() > p2.first->GetName()) {
44  result = true;
45  }
46  }
47  }
48  return result;
49  }
50 
51  bool sortByName(const std::pair<G4LogicalVolume*, const dd4hep::SpecPar*>& p1,
52  const std::pair<G4LogicalVolume*, const dd4hep::SpecPar*>& p2) {
53  bool result = false;
54  if (p1.first->GetName() > p2.first->GetName()) {
55  result = true;
56  }
57  return result;
58  }
59 } // namespace
60 
62  : map_(map), keywordRegion_("CMSCutsRegion"), verbosity_(verb), protonCut_(pcut) {
63  initialize();
64 }
65 
66 DDG4ProductionCuts::DDG4ProductionCuts(const dd4hep::SpecParRegistry* specPars,
67  const dd4hep::sim::Geant4GeometryMaps::VolumeMap* map,
68  int verb,
69  bool pcut)
70  : dd4hepMap_(map), specPars_(specPars), keywordRegion_("CMSCutsRegion"), verbosity_(verb), protonCut_(pcut) {
72 }
73 
76  // sort all root volumes - to get the same sequence at every run of the application.
77  // (otherwise, the sequence will depend on the pointer (memory address) of the
78  // involved objects, because 'new' does no guarantee that you allways get a
79  // higher (or lower) address when allocating an object of the same type ...
80  sort(vec_.begin(), vec_.end(), &dd_is_greater);
81  if (verbosity_ > 0) {
82  edm::LogVerbatim("Geometry") << " DDG4ProductionCuts : got " << vec_.size() << " region roots.\n"
83  << " DDG4ProductionCuts : List of all roots:";
84  for (auto const& vv : vec_)
85  edm::LogVerbatim("Geometry") << " " << vv.first->GetName() << " : " << vv.second.name();
86  }
87 
88  // Now generate all the regions
89  std::string curName = "";
90  std::string regionName = "";
91  G4Region* region = nullptr;
92  G4RegionStore* store = G4RegionStore::GetInstance();
93  for (auto const& vv : vec_) {
94  unsigned int num = map_->toString(keywordRegion_, vv.second, regionName);
95  edm::LogVerbatim("Geometry") << " num " << num << " regionName: " << regionName << ", the store of size "
96  << store->size();
97  if (num != 1) {
98  throw cms::Exception("SimG4CoreGeometry", " DDG4ProductionCuts::initialize: Problem with Region tags.");
99  }
100  if (regionName != curName) {
101  edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : regionName " << regionName << ", the store of size "
102  << store->size();
103  region = store->FindOrCreateRegion(regionName);
104  edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : region " << region->GetName();
105  if (!region) {
106  throw cms::Exception("SimG4CoreGeometry", " DDG4ProductionCuts::initialize: Problem with Region tags.");
107  }
108  curName = regionName;
109  edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : new G4Region " << vv.first->GetName();
110  setProdCuts(vv.second, region);
111  }
112 
113  region->AddRootLogicalVolume(vv.first);
114 
115  if (verbosity_ > 0)
116  edm::LogVerbatim("Geometry") << " added " << vv.first->GetName() << " to region " << region->GetName();
117  }
118 }
119 
121  dd4hep::SpecParRefs specs;
122  specPars_->filter(specs, keywordRegion_);
123 
124  for (auto const& it : *dd4hepMap_) {
125  for (auto const& fit : specs) {
126  for (auto const& pit : fit.second->paths) {
127  const std::string_view selection = dd4hep::dd::realTopName(pit);
128  const std::string_view name = dd4hep::dd::noNamespace(it.first.name());
129  if (!(dd4hep::dd::isRegex(selection))
130  ? dd4hep::dd::compareEqual(name, selection)
131  : std::regex_match(name.begin(), name.end(), std::regex(std::string(selection)))) {
132  dd4hepVec_.emplace_back(std::make_pair<G4LogicalVolume*, const dd4hep::SpecPar*>(&*it.second, &*fit.second));
133  }
134  }
135  }
136  }
137  // sort all root volumes - to get the same sequence at every run of the application.
138  sort(begin(dd4hepVec_), end(dd4hepVec_), &sortByName);
139 
140  // Now generate all the regions
141  for (auto const& it : dd4hepVec_) {
142  auto regName = it.second->strValue(keywordRegion_);
143  G4Region* region = G4RegionStore::GetInstance()->FindOrCreateRegion({regName.data(), regName.size()});
144  region->AddRootLogicalVolume(it.first);
145  edm::LogVerbatim("Geometry") << it.first->GetName() << ": " << it.second->strValue(keywordRegion_);
146  edm::LogVerbatim("Geometry") << " MakeRegions: added " << it.first->GetName() << " to region " << region->GetName();
147  edm::LogVerbatim("Geometry").log([&](auto& log) {
148  for (auto const& sit : it.second->spars) {
149  log << sit.first << " = " << sit.second[0] << "\n";
150  }
151  });
152  setProdCuts(it.second, region);
153  }
154 
155  if (verbosity_ > 0) {
156  edm::LogVerbatim("SimG4CoreGeometry") << " DDG4ProductionCuts (New) : starting\n"
157  << " DDG4ProductionCuts : Got " << dd4hepVec_.size() << " region roots.\n"
158  << " DDG4ProductionCuts : List of all roots:";
159  for (size_t jj = 0; jj < dd4hepVec_.size(); ++jj)
160  edm::LogVerbatim("SimG4CoreGeometry") << " DDG4ProductionCuts : root=" << dd4hepVec_[jj].first->GetName()
161  << " , " << dd4hepVec_[jj].second->paths.at(0);
162  }
163 }
164 
166  //
167  // search for production cuts
168  // you must have four of them: e+ e- gamma proton
169  //
170  double gammacut = 0.0;
171  double electroncut = 0.0;
172  double positroncut = 0.0;
173  double protoncut = 0.0;
174  int temp = map_->toDouble("ProdCutsForGamma", lpart, gammacut);
175  if (temp != 1) {
176  throw cms::Exception(
177  "SimG4CorePhysics",
178  " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForGamma.");
179  }
180  temp = map_->toDouble("ProdCutsForElectrons", lpart, electroncut);
181  if (temp != 1) {
182  throw cms::Exception(
183  "SimG4CorePhysics",
184  " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForElectrons.");
185  }
186  temp = map_->toDouble("ProdCutsForPositrons", lpart, positroncut);
187  if (temp != 1) {
188  throw cms::Exception(
189  "SimG4CorePhysics",
190  " DDG4ProductionCuts::setProdCuts: Problem with Region tags - no/more than one ProdCutsForPositrons.");
191  }
192  temp = map_->toDouble("ProdCutsForProtons", lpart, protoncut);
193  if (temp == 0) {
194  // There is no ProdCutsForProtons set in XML,
195  // check if it's a legacy geometry scenario without it
196  if (protonCut_) {
197  protoncut = electroncut;
198  } else {
199  protoncut = 0.;
200  }
201  } else if (temp != 1) {
202  throw cms::Exception(
203  "SimG4CorePhysics",
204  " DDG4ProductionCuts::setProdCuts: Problem with Region tags - more than one ProdCutsForProtons.");
205  }
206 
207  //
208  // Create and fill production cuts
209  //
210  G4ProductionCuts* prodCuts = region->GetProductionCuts();
211  if (!prodCuts) {
212  prodCuts = new G4ProductionCuts();
213  region->SetProductionCuts(prodCuts);
214  }
215  prodCuts->SetProductionCut(gammacut, idxG4GammaCut);
216  prodCuts->SetProductionCut(electroncut, idxG4ElectronCut);
217  prodCuts->SetProductionCut(positroncut, idxG4PositronCut);
218  prodCuts->SetProductionCut(protoncut, idxG4ProtonCut);
219  if (verbosity_ > 0) {
220  edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : Setting cuts for " << region->GetName()
221  << "\n Electrons: " << electroncut << "\n Positrons: " << positroncut
222  << "\n Gamma : " << gammacut << "\n Proton : " << protoncut;
223  }
224 }
225 
226 void DDG4ProductionCuts::setProdCuts(const dd4hep::SpecPar* spec, G4Region* region) {
227  //
228  // Create and fill production cuts
229  //
230  G4ProductionCuts* prodCuts = region->GetProductionCuts();
231  if (!prodCuts) {
232  //
233  // search for production cuts
234  // you must have four of them: e+ e- gamma proton
235  //
236  double gammacut = spec->dblValue("ProdCutsForGamma") / MM_2_CM;
237  double electroncut = spec->dblValue("ProdCutsForElectrons") / MM_2_CM;
238  double positroncut = spec->dblValue("ProdCutsForPositrons") / MM_2_CM;
239  double protoncut = spec->dblValue("ProdCutsForProtons") / MM_2_CM;
240  if (protoncut == 0) {
241  protoncut = electroncut;
242  }
243 
244  prodCuts = new G4ProductionCuts();
245  region->SetProductionCuts(prodCuts);
246 
247  prodCuts->SetProductionCut(gammacut, idxG4GammaCut);
248  prodCuts->SetProductionCut(electroncut, idxG4ElectronCut);
249  prodCuts->SetProductionCut(positroncut, idxG4PositronCut);
250  prodCuts->SetProductionCut(protoncut, idxG4ProtonCut);
251  if (verbosity_ > 0) {
252  edm::LogVerbatim("Geometry") << "DDG4ProductionCuts : Setting cuts for " << region->GetName()
253  << "\n Electrons: " << electroncut << "\n Positrons: " << positroncut
254  << "\n Gamma : " << gammacut << "\n Proton : " << protoncut;
255  }
256  } else {
257  if (verbosity_ > 0) {
258  edm::LogVerbatim("Geometry")
259  << "DDG4ProductionCuts : Cuts are already set for " << region->GetName()
260  << "\n Electrons: " << region->GetProductionCuts()->GetProductionCut(idxG4ElectronCut)
261  << "\n Positrons: " << region->GetProductionCuts()->GetProductionCut(idxG4PositronCut)
262  << "\n Gamma : " << region->GetProductionCuts()->GetProductionCut(idxG4GammaCut)
263  << "\n Proton : " << region->GetProductionCuts()->GetProductionCut(idxG4ProtonCut);
264  }
265  }
266 }
DDMapper::toDouble
unsigned int toDouble(const std::string &name, const KeyType &key, double &value, unsigned int pos=0) const
returns the number specific parameters named 'name' and the corrsponding double
Definition: DDMapper.h:111
DDG4ProductionCuts::verbosity_
const int verbosity_
Definition: DDG4ProductionCuts.h:44
MessageLogger.h
DDMapper< G4LogicalVolume *, DDLogicalPart >
DDG4ProductionCuts::DDG4ProductionCuts
DDG4ProductionCuts(const G4LogicalVolumeToDDLogicalPartMap *, int, bool)
Definition: DDG4ProductionCuts.cc:61
HistogramManager_cfi.specs
specs
Definition: HistogramManager_cfi.py:80
DDG4ProductionCuts::dd4hepMap_
const dd4hep::sim::Geant4GeometryMaps::VolumeMap * dd4hepMap_
Definition: DDG4ProductionCuts.h:37
groupFilesInBlocks.temp
list temp
Definition: groupFilesInBlocks.py:142
DDG4ProductionCuts::protonCut_
const bool protonCut_
Definition: DDG4ProductionCuts.h:45
DDG4ProductionCuts::map_
const G4LogicalVolumeToDDLogicalPartMap * map_
Definition: DDG4ProductionCuts.h:28
DDG4ProductionCuts.h
mps_fire.end
end
Definition: mps_fire.py:242
corrVsCorr.selection
selection
main part
Definition: corrVsCorr.py:100
p2
double p2[4]
Definition: TauolaWrapper.h:90
AlCaHLTBitMon_QueryRunRegistry.string
string
Definition: AlCaHLTBitMon_QueryRunRegistry.py:256
geant_units::operators::convertCmToMm
constexpr NumType convertCmToMm(NumType centimeters)
Definition: GeantUnits.h:68
DDG4ProductionCuts::keywordRegion_
const std::string keywordRegion_
Definition: DDG4ProductionCuts.h:43
GeantUnits.h
DDLogicalPart
A DDLogicalPart aggregates information concerning material, solid and sensitveness ....
Definition: DDLogicalPart.h:93
MM_2_CM
#define MM_2_CM
Definition: DDG4ProductionCuts.cc:23
HLT_FULL_cff.region
region
Definition: HLT_FULL_cff.py:84949
DDG4ProductionCuts::dd4hepVec_
std::vector< std::pair< G4LogicalVolume *, const dd4hep::SpecPar * > > dd4hepVec_
Definition: DDG4ProductionCuts.h:38
DDG4ProductionCuts::specPars_
const dd4hep::SpecParRegistry * specPars_
Definition: DDG4ProductionCuts.h:39
p1
double p1[4]
Definition: TauolaWrapper.h:89
DDMapper::toString
unsigned int toString(const std::string &name, const KeyType &key, std::string &value, unsigned int pos=0) const
same as toDouble but for std::string-valued values of named parameters
Definition: DDMapper.h:172
DDLogicalPart.h
EgammaValidation_cff.num
num
Definition: EgammaValidation_cff.py:34
DDG4ProductionCuts::vec_
G4LogicalVolumeToDDLogicalPartMap::Vector vec_
Definition: DDG4ProductionCuts.h:29
edm::LogVerbatim
Log< level::Info, true > LogVerbatim
Definition: MessageLogger.h:128
Exception
Definition: hltDiff.cc:246
findQualityFiles.jj
string jj
Definition: findQualityFiles.py:188
DDG4ProductionCuts::initialize
void initialize()
Definition: DDG4ProductionCuts.cc:74
DDMapper::all
Vector all(const std::string &name, const std::string &value) const
get all std::mapped instances which have a specific 'name' with value 'value'
Definition: DDMapper.h:199
Skims_PA_cff.name
name
Definition: Skims_PA_cff.py:17
Exception.h
dqm-mbProfile.log
log
Definition: dqm-mbProfile.py:17
mps_fire.result
result
Definition: mps_fire.py:311
genParticles_cff.map
map
Definition: genParticles_cff.py:11
DDG4ProductionCuts::setProdCuts
void setProdCuts(const DDLogicalPart, G4Region *)
Definition: DDG4ProductionCuts.cc:165
edm::Log
Definition: MessageLogger.h:70
DDG4ProductionCuts::dd4hepInitialize
void dd4hepInitialize()
Definition: DDG4ProductionCuts.cc:120
fit
Definition: CombinedChiSquaredLikelihood.h:6