SimplifiedGeometryPropagator/src/SimplifiedGeometryFactory.cc

0001 #include "FastSimulation/SimplifiedGeometryPropagator/interface/SimplifiedGeometryFactory.h"
0002
0003 #include "FWCore/Utilities/interface/Exception.h"
0004 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0005 #include "FastSimulation/SimplifiedGeometryPropagator/interface/BarrelSimplifiedGeometry.h"
0006 #include "FastSimulation/SimplifiedGeometryPropagator/interface/ForwardSimplifiedGeometry.h"
0007 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0008 #include "MagneticField/Engine/interface/MagneticField.h"
0009
0010 #include <TH1F.h>
0011 #include <cctype>
0012 #include <memory>
0013
0014 fastsim::SimplifiedGeometryFactory::SimplifiedGeometryFactory(
0015     const GeometricSearchTracker *geometricSearchTracker,
0016     const MagneticField &magneticField,
0017     const std::map<std::string, fastsim::InteractionModel *> &interactionModelMap,
0018     double magneticFieldHistMaxR,
0019     double magneticFieldHistMaxZ)
0020     : geometricSearchTracker_(geometricSearchTracker),
0021       magneticField_(&magneticField),
0022       interactionModelMap_(&interactionModelMap),
0023       magneticFieldHistMaxR_(magneticFieldHistMaxR),
0024       magneticFieldHistMaxZ_(magneticFieldHistMaxZ) {
0025   // naming convention for barrel DetLayer lists
0026   barrelDetLayersMap_["BPix"] = &geometricSearchTracker_->pixelBarrelLayers();
0027   barrelDetLayersMap_["TIB"] = &geometricSearchTracker_->tibLayers();
0028   barrelDetLayersMap_["TOB"] = &geometricSearchTracker_->tobLayers();
0029
0030   // naming convention for forward DetLayer lists
0031   forwardDetLayersMap_["negFPix"] = &geometricSearchTracker_->negPixelForwardLayers();
0032   forwardDetLayersMap_["posFPix"] = &geometricSearchTracker_->posPixelForwardLayers();
0033   forwardDetLayersMap_["negTID"] = &geometricSearchTracker_->negTidLayers();
0034   forwardDetLayersMap_["posTID"] = &geometricSearchTracker_->posTidLayers();
0035   forwardDetLayersMap_["negTEC"] = &geometricSearchTracker_->negTecLayers();
0036   forwardDetLayersMap_["posTEC"] = &geometricSearchTracker_->posTecLayers();
0037 }
0038
0039 std::unique_ptr<fastsim::BarrelSimplifiedGeometry> fastsim::SimplifiedGeometryFactory::createBarrelSimplifiedGeometry(
0040     const edm::ParameterSet &cfg) const {
0041   std::unique_ptr<fastsim::SimplifiedGeometry> layer = createSimplifiedGeometry(BARREL, cfg);
0042   return std::unique_ptr<fastsim::BarrelSimplifiedGeometry>(
0043       static_cast<fastsim::BarrelSimplifiedGeometry *>(layer.release()));
0044 }
0045
0046 std::unique_ptr<fastsim::ForwardSimplifiedGeometry> fastsim::SimplifiedGeometryFactory::createForwardSimplifiedGeometry(
0047     const fastsim::SimplifiedGeometryFactory::LayerType layerType, const edm::ParameterSet &cfg) const {
0048   if (layerType != NEGFWD && layerType != POSFWD) {
0049     throw cms::Exception("fastsim::SimplifiedGeometry::createForwardLayer")
0050         << " called with forbidden layerType. Allowed layerTypes are NEGFWD and POSFWD";
0051   }
0052   std::unique_ptr<fastsim::SimplifiedGeometry> layer = createSimplifiedGeometry(layerType, cfg);
0053   return std::unique_ptr<fastsim::ForwardSimplifiedGeometry>(
0054       static_cast<fastsim::ForwardSimplifiedGeometry *>(layer.release()));
0055 }
0056
0057 std::unique_ptr<fastsim::SimplifiedGeometry> fastsim::SimplifiedGeometryFactory::createSimplifiedGeometry(
0058     const fastsim::SimplifiedGeometryFactory::LayerType layerType, const edm::ParameterSet &cfg) const {
0059   // some flags for internal usage
0060   bool isForward = true;
0061   bool isOnPositiveSide = false;
0062   if (layerType == BARREL) {
0063     isForward = false;
0064   } else if (layerType == POSFWD) {
0065     isOnPositiveSide = true;
0066   }
0067
0068   // -------------------------------
0069   // extract DetLayer (i.e. full geometry of tracker modules)
0070   // -------------------------------
0071
0072   std::string detLayerName = cfg.getUntrackedParameter<std::string>("activeLayer", "");
0073   const DetLayer *detLayer = nullptr;
0074
0075   if (!detLayerName.empty() && geometricSearchTracker_) {
0076     if (isForward) {
0077       detLayerName = (isOnPositiveSide ? "pos" : "neg") + detLayerName;
0078     }
0079     detLayer = getDetLayer(detLayerName, *geometricSearchTracker_);
0080   }
0081
0082   // ------------------------------
0083   // radius / z of layers
0084   // ------------------------------
0085
0086   // first try to get it from the configuration
0087   double position = 0;
0088   std::string positionParameterName = (isForward ? "z" : "radius");
0089   if (cfg.exists(positionParameterName)) {
0090     position = fabs(cfg.getUntrackedParameter<double>(positionParameterName));
0091     if (isForward && !isOnPositiveSide) {
0092       position = -position;
0093     }
0094   }
0095   // then try extracting from detLayer
0096   else if (detLayer) {
0097     if (isForward) {
0098       position = static_cast<ForwardDetLayer const *>(detLayer)->surface().position().z();
0099     } else {
0100       position = static_cast<BarrelDetLayer const *>(detLayer)->specificSurface().radius();
0101     }
0102   }
0103   // then throw error
0104   else {
0105     std::string cfgString;
0106     cfg.allToString(cfgString);
0107     throw cms::Exception("fastsim::SimplifiedGeometry")
0108         << "Cannot extract a " << (isForward ? "position" : "radius") << " for this "
0109         << (isForward ? "forward" : "barrel") << " layer:\n"
0110         << cfgString;
0111   }
0112
0113   // -----------------------------
0114   // create the layers
0115   // -----------------------------
0116
0117   std::unique_ptr<fastsim::SimplifiedGeometry> layer;
0118   if (isForward) {
0119     layer = std::make_unique<fastsim::ForwardSimplifiedGeometry>(position);
0120   } else {
0121     layer = std::make_unique<fastsim::BarrelSimplifiedGeometry>(position);
0122   }
0123   layer->detLayer_ = detLayer;
0124
0125   // -----------------------------
0126   // thickness histogram
0127   // -----------------------------
0128
0129   // Get limits
0130   const std::vector<double> &limits = cfg.getUntrackedParameter<std::vector<double> >("limits");
0131   // and check order.
0132   for (unsigned index = 1; index < limits.size(); index++) {
0133     if (limits[index] < limits[index - 1]) {
0134       std::string cfgString;
0135       cfg.allToString(cfgString);
0136       throw cms::Exception("fastsim::SimplifiedGeometryFactory")
0137           << "limits must be provided in increasing order. error in:\n"
0138           << cfgString;
0139     }
0140   }
0141   // Get thickness values
0142   const std::vector<double> &thickness = cfg.getUntrackedParameter<std::vector<double> >("thickness");
0143   // and check compatibility with limits
0144   if (limits.size() < 2 || thickness.size() != limits.size() - 1) {
0145     std::string cfgString;
0146     cfg.allToString(cfgString);
0147     throw cms::Exception("fastim::SimplifiedGeometryFactory")
0148         << "layer thickness and limits not configured properly! error in:" << cfgString;
0149   }
0150   // create the histogram
0151   layer->thicknessHist_ = std::make_unique<TH1F>("h", "h", limits.size() - 1, &limits[0]);
0152   layer->thicknessHist_->SetDirectory(nullptr);
0153   for (unsigned i = 1; i < limits.size(); ++i) {
0154     layer->thicknessHist_->SetBinContent(i, thickness[i - 1]);
0155   }
0156
0157   // -----------------------------
0158   // nuclear interaction thickness factor
0159   // -----------------------------
0160
0161   layer->nuclearInteractionThicknessFactor_ =
0162       cfg.getUntrackedParameter<double>("nuclearInteractionThicknessFactor", 1.);
0163
0164   // -----------------------------
0165   // magnetic field
0166   // -----------------------------
0167
0168   layer->magneticFieldHist_ =
0169       std::make_unique<TH1F>("h", "h", 100, 0., isForward ? magneticFieldHistMaxR_ : magneticFieldHistMaxZ_);
0170   layer->magneticFieldHist_->SetDirectory(nullptr);
0171   for (int i = 1; i <= 101; i++) {
0172     GlobalPoint point = isForward ? GlobalPoint(layer->magneticFieldHist_->GetXaxis()->GetBinCenter(i), 0., position)
0173                                   : GlobalPoint(position, 0., layer->magneticFieldHist_->GetXaxis()->GetBinCenter(i));
0174     layer->magneticFieldHist_->SetBinContent(i, magneticField_->inTesla(point).z());
0175   }
0176
0177   // -----------------------------
0178   // list of interaction models
0179   // -----------------------------
0180
0181   std::vector<std::string> interactionModelLabels =
0182       cfg.getUntrackedParameter<std::vector<std::string> >("interactionModels");
0183   for (const auto &label : interactionModelLabels) {
0184     std::map<std::string, fastsim::InteractionModel *>::const_iterator interactionModel =
0185         interactionModelMap_->find(label);
0186     if (interactionModel == interactionModelMap_->end()) {
0187       throw cms::Exception("fastsim::SimplifiedGeometryFactory") << "unknown interaction model '" << label << "'";
0188     }
0189     layer->interactionModels_.push_back(interactionModel->second);
0190   }
0191
0192   // -----------------------------
0193   // Hack to interface "old" calorimetry with "new" propagation in tracker
0194   // -----------------------------
0195
0196   if (cfg.exists("caloType")) {
0197     std::string caloType = cfg.getUntrackedParameter<std::string>("caloType");
0198
0199     if (caloType == "PRESHOWER1") {
0200       layer->setCaloType(SimplifiedGeometry::PRESHOWER1);
0201     } else if (caloType == "PRESHOWER2") {
0202       layer->setCaloType(SimplifiedGeometry::PRESHOWER2);
0203     } else if (caloType == "ECAL") {
0204       layer->setCaloType(SimplifiedGeometry::ECAL);
0205     } else if (caloType == "HCAL") {
0206       layer->setCaloType(SimplifiedGeometry::HCAL);
0207     } else if (caloType == "VFCAL") {
0208       layer->setCaloType(SimplifiedGeometry::VFCAL);
0209     } else {
0210       throw cms::Exception("fastsim::SimplifiedGeometryFactory")
0211           << "unknown caloType '" << caloType << "' (defined PRESHOWER1, PRESHOWER2, ECAL, HCAL, VFCAL)";
0212     }
0213   }
0214
0215   // -----------------------------
0216   // and return the layer!
0217   // -----------------------------
0218
0219   return layer;
0220 }
0221
0222 const DetLayer *fastsim::SimplifiedGeometryFactory::getDetLayer(
0223     const std::string &detLayerName, const GeometricSearchTracker &geometricSearchTracker) const {
0224   if (detLayerName.empty()) {
0225     return nullptr;
0226   }
0227
0228   // obtain the index from the detLayerName
0229   unsigned pos = detLayerName.size();
0230   while (isdigit(detLayerName[pos - 1])) {
0231     pos -= 1;
0232   }
0233   if (pos == detLayerName.size()) {
0234     throw cms::Exception("fastsim::SimplifiedGeometry::getDetLayer")
0235         << "last part of detLayerName must be index of DetLayer in list. Error in detLayerName" << detLayerName
0236         << std::endl;
0237   }
0238   int index = atoi(detLayerName.substr(pos).c_str());
0239   std::string detLayerListName = detLayerName.substr(0, pos);
0240
0241   try {
0242     // try to find the detLayer in the barrel map
0243     if (barrelDetLayersMap_.find(detLayerListName) != barrelDetLayersMap_.end()) {
0244       auto detLayerList = barrelDetLayersMap_.find(detLayerListName)->second;
0245       return detLayerList->at(index - 1);  // use at, to provoce the throwing of an error in case of a bad index
0246     }
0247
0248     // try to find the detLayer in the forward map
0249     else if (forwardDetLayersMap_.find(detLayerListName) != forwardDetLayersMap_.end()) {
0250       auto detLayerList = forwardDetLayersMap_.find(detLayerListName)->second;
0251       return detLayerList->at(index - 1);  // use at, to provoce the throwing of an error in case of a bad index
0252     }
0253     // throw an error
0254     else {
0255       throw cms::Exception("fastsim::SimplifiedGeometry::getDetLayer")
0256           << " could not find list of detLayers corresponding to detLayerName " << detLayerName << std::endl;
0257     }
0258   } catch (const std::out_of_range &error) {
0259     throw cms::Exception("fastsim::SimplifiedGeometry::getDetLayer")
0260         << " index out of range for detLayerName: " << detLayerName << " " << error.what() << std::endl;
0261   }
0262 }