From af2bb20d166ba274a4d2cf888fd2d12c082bbf5a Mon Sep 17 00:00:00 2001 From: Kevin Pedro Date: Wed, 26 Sep 2018 10:02:18 -0500 Subject: [PATCH 1/2] squashed #23799 (New pixel RecHit infrastructure for FastSim (for Phases 0, 1, 2, and beyond)) --- .../interface/PixelResolutionHistograms.h | 135 ++++ .../interface/PixelTemplateSmearerBase.h | 57 +- .../interface/TrackingRecHitAlgorithm.h | 12 +- .../PixelBarrelTemplateSmearerPlugin.cc | 130 ---- .../PixelForwardTemplateSmearerPlugin.cc | 150 ----- .../plugins/PixelTemplateSmearerPlugin.cc | 51 ++ .../plugins/TrackingRecHitProducer.cc | 33 +- .../python/PixelPluginsPhase0_cfi.py | 59 ++ .../python/PixelPluginsPhase1_cfi.py | 191 ++++++ .../python/PixelPluginsPhase2_cfi.py | 46 ++ .../python/TrackingRecHitProducer_cfi.py | 66 +- .../src/PixelResolutionHistograms.cc | 580 ++++++++++++++++++ .../src/PixelTemplateSmearerBase.cc | 534 +++++++++------- .../src/TrackingRecHitAlgorithm.cc | 9 + 14 files changed, 1485 insertions(+), 568 deletions(-) create mode 100644 FastSimulation/TrackingRecHitProducer/interface/PixelResolutionHistograms.h delete mode 100644 FastSimulation/TrackingRecHitProducer/plugins/PixelBarrelTemplateSmearerPlugin.cc delete mode 100644 FastSimulation/TrackingRecHitProducer/plugins/PixelForwardTemplateSmearerPlugin.cc create mode 100644 FastSimulation/TrackingRecHitProducer/plugins/PixelTemplateSmearerPlugin.cc create mode 100644 FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase0_cfi.py create mode 100644 FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase1_cfi.py create mode 100644 FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase2_cfi.py create mode 100644 FastSimulation/TrackingRecHitProducer/src/PixelResolutionHistograms.cc diff --git a/FastSimulation/TrackingRecHitProducer/interface/PixelResolutionHistograms.h b/FastSimulation/TrackingRecHitProducer/interface/PixelResolutionHistograms.h new file mode 100644 index 0000000000000..9aeabb3c70424 --- /dev/null +++ b/FastSimulation/TrackingRecHitProducer/interface/PixelResolutionHistograms.h @@ -0,0 +1,135 @@ +#ifndef FastSimulation_TrackingRecHitProducer_PixelResolutionHistograms_h +#define FastSimulation_TrackingRecHitProducer_PixelResolutionHistograms_h 1 + +class TFile; +class TH1F; +class TH2F; +class TAxis; +class RandomEngineAndDistribution; +class SimpleHistogramGenerator; + +#include + +/// #define COTBETA_HIST_MAX 30 +/// #define COTALPHA_HIST_MAX 20 +/// #define QBIN_HIST_MAX 4 + +static constexpr unsigned int COTBETA_HIST_MAX = 30; +static constexpr unsigned int COTALPHA_HIST_MAX = 20; +static constexpr unsigned int QBIN_HIST_MAX = 4; + +class PixelResolutionHistograms { +public: + + //--- Constructor to use when generating resolution histograms. + // We make empty histograms (which we own), but generator pointers + // remain null. + // + PixelResolutionHistograms( std::string filename, // ROOT file for histograms + std::string rootdir, // Subdirectory in the file, "" if none + std::string descTitle, // Descriptive title + unsigned int detType, // Where we are... (&&& do we need this?) + double cotbetaBinWidth, // cot(beta) : bin width, + double cotbetaLowEdge, // : low endpoint, + int cotbetaBins, // : # of bins + double cotalphaBinWidth, // cot(alpha): bin width, + double cotalphaLowEdge, // : low endpoint, + int cotalphaBins ); // : # of bins + //int qbinWidth, + //int qbins ) + + + //--- Constructor to use when reading the histograms from a file (e.g. when + // inside a running FastSim job). We get the histograms from a + // ROOT file, and we do *not* own them. But we do own the + // generators. + // + PixelResolutionHistograms( std::string filename, // ROOT file for histograms + std::string rootdir = "", // ROOT dir, "" if none + int detType = -1, // default: read from ROOT file. + bool ignore_multi = false, // Forward Big is always single + bool ignore_single = false, // Edge does not need single pixels + bool ignore_qBin = false ); // qBin histograms not used right now (future expansion) + + + //--- Destructor (virtual, just in case) + virtual ~PixelResolutionHistograms(); + + //--- Status after construction (esp.loading from file). Non-zero if there + // were problems. + inline int status() { return status_ ; } + + //--- Fill one entry in one resolution histogram. Use when making histograms. + int Fill( double dx, double dy, // the difference wrt true hit + double cotalpha, double cotbeta, // cotangent of local angles + int qbin, // Qbin = category for how much charge we have + int nxpix, int nypix ); // length of cluster along x,y (only care if ==1 or not) + + + //--- Get generators, for resolution in X and Y. Use in FastSim. + const SimpleHistogramGenerator * getGeneratorX( double cotalpha, + double cotbeta, + int qbin, + bool singlex ); + + const SimpleHistogramGenerator * getGeneratorY( double cotalpha, + double cotbeta, + int qbin, + bool singley ); + + + private: + // Do we own the histograms, or not? + bool weOwnHistograms_ ; + + // Where we are. + unsigned int detType_ ; // 1 for barrel, 0 for forward /// May not need this? + + // Resolution binning + double cotbetaBinWidth_ ; + double cotbetaLowEdge_ ; + int cotbetaBins_ ; + double cotalphaBinWidth_ ; + double cotalphaLowEdge_ ; + int cotalphaBins_ ; + int qbinWidth_ ; + int qbins_ ; + + // The dummy histogram to hold the binning, and the two cached axes. + TH2F * binningHisto_ ; + TAxis * cotbetaAxis_ ; + TAxis * cotalphaAxis_ ; + + // Resolution histograms. I (Petar) tried to dynamically allocate + // these histograms, but all possible implementations were somewhat + // complicated, which would make the code harder to understand, + // debug, and thus support in the long term. Since we are here only + // booking pointers of histograms, we will instead book larger + // matrices, and leave them partially empty. But who cares -- the + // wasted memory of a few hundred null pointers is negligible. + // + // The constructor will then fill only the first cotbetaBins_ x + // cotalphaBins_ x qbinBins_ histograms in the matrix, and we'll + // ignore the rest. + // + TH1F * resMultiPixelXHist_ [ COTBETA_HIST_MAX ][ COTALPHA_HIST_MAX ][ QBIN_HIST_MAX ] ; + TH1F * resSinglePixelXHist_ [ COTBETA_HIST_MAX ][ COTALPHA_HIST_MAX ]; + TH1F * resMultiPixelYHist_ [ COTBETA_HIST_MAX ][ COTALPHA_HIST_MAX ][ QBIN_HIST_MAX ]; + TH1F * resSinglePixelYHist_ [ COTBETA_HIST_MAX ][ COTALPHA_HIST_MAX ]; + TH1F * qbinHist_ [ COTBETA_HIST_MAX ][ COTALPHA_HIST_MAX ]; + + // File with histograms to load. + std::unique_ptr file_ ; + + // Status of loading. Check if there were errors. + int status_ ; + + // Identical binning and parameterization for FastSim generators. + SimpleHistogramGenerator * resMultiPixelXGen_ [ COTBETA_HIST_MAX ][ COTALPHA_HIST_MAX ][ QBIN_HIST_MAX ] ; + SimpleHistogramGenerator * resSinglePixelXGen_ [ COTBETA_HIST_MAX ][ COTALPHA_HIST_MAX ]; + SimpleHistogramGenerator * resMultiPixelYGen_ [ COTBETA_HIST_MAX ][ COTALPHA_HIST_MAX ][ QBIN_HIST_MAX ]; + SimpleHistogramGenerator * resSinglePixelYGen_ [ COTBETA_HIST_MAX ][ COTALPHA_HIST_MAX ]; + SimpleHistogramGenerator * qbinGen_ [ COTBETA_HIST_MAX ][ COTALPHA_HIST_MAX ]; + +}; +#endif diff --git a/FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h b/FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h index 0aa7f7b3e987c..387990ad00113 100644 --- a/FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h +++ b/FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h @@ -27,7 +27,7 @@ #include "DataFormats/GeometryVector/interface/Point3DBase.h" #include "DataFormats/GeometrySurface/interface/LocalError.h" -// STL +// STL. needed for uniq_ptr<> #include #include #include @@ -35,6 +35,7 @@ class TFile; class RandomEngineAndDistribution; class SimpleHistogramGenerator; +class PixelResolutionHistograms; class PixelTemplateSmearerBase: public TrackingRecHitAlgorithm @@ -47,48 +48,54 @@ class PixelTemplateSmearerBase: }; protected: - bool mergeHitsOn; - std::vector< SiPixelTemplateStore > thePixelTemp_; - int templateId; - - bool isFlipped(const PixelGeomDetUnit* theDet) const; - //isForward, true for forward, false for barrel - bool isForward; + bool mergeHitsOn = false; // if true then see if neighboring hits might merge + + //--- Template DB Object(s) + const SiPixelTemplateDBObject * pixelTemplateDBObject_ = nullptr; // needed for template<-->DetId map. + std::vector< SiPixelTemplateStore > thePixelTemp_ ; // our own template storage + std::vector< SiPixelTemplateStore > & thePixelTempRef = thePixelTemp_; // points to the one we will use. + int templateId = -1; - double rescotAlpha_binMin , rescotAlpha_binWidth; - unsigned int rescotAlpha_binN; - double rescotBeta_binMin , rescotBeta_binWidth; - unsigned int rescotBeta_binN; - int resqbin_binMin, resqbin_binWidth; - unsigned int resqbin_binN; + //--- Flag to tell us whether we are in barrel or in forward. + // This is needed since the parameterization is slightly + // different for forward, since all forward detectors cover + // a smaller range of local incidence angles and thus + // the clusters are shorter and have less charge. + bool isBarrel; - - std::map theXHistos; - std::map theYHistos; - - std::unique_ptr theEdgePixelResolutionFile; + //--- The histogram storage containers. + std::shared_ptr theEdgePixelResolutions; std::string theEdgePixelResolutionFileName; - std::unique_ptr theBigPixelResolutionFile; + + std::shared_ptr theBigPixelResolutions; std::string theBigPixelResolutionFileName; - std::unique_ptr theRegularPixelResolutionFile; + + std::shared_ptr theRegularPixelResolutions; std::string theRegularPixelResolutionFileName; + + //--- Files with hit merging information: std::unique_ptr theMergingProbabilityFile; std::string theMergingProbabilityFileName; + std::unique_ptr theMergedPixelResolutionXFile; std::string theMergedPixelResolutionXFileName; - std::unique_ptr theMergedPixelResolutionYFile; - std::string theMergedPixelResolutionYFileName; - unsigned int theLayer; + std::unique_ptr theMergedPixelResolutionYFile; + std::string theMergedPixelResolutionYFileName; + public: - explicit PixelTemplateSmearerBase( const std::string& name, const edm::ParameterSet& config, edm::ConsumesCollector& consumesCollector ); ~PixelTemplateSmearerBase() override; TrackingRecHitProductPtr process(TrackingRecHitProductPtr product) const override; + // void beginEvent(edm::Event& event, const edm::EventSetup& eventSetup) override; + void beginRun(edm::Run const& run, const edm::EventSetup& eventSetup, + const SiPixelTemplateDBObject * pixelTemplateDBObjectPtr, + std::vector< SiPixelTemplateStore > & tempStoreRef ) override; + // void endEvent(edm::Event& event, const edm::EventSetup& eventSetup) override; //--- Process all unmerged hits. Calls smearHit() for each. TrackingRecHitProductPtr processUnmergedHits( diff --git a/FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithm.h b/FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithm.h index edb81db3d0109..9c381028c73fe 100644 --- a/FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithm.h +++ b/FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithm.h @@ -3,11 +3,15 @@ #include "DataFormats/TrackerCommon/interface/TrackerTopology.h" #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" - +#include "FWCore/Framework/interface/ProducerBase.h" #include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitProduct.h" #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h" +// Pixel-related stuff: +#include "CondFormats/SiPixelObjects/interface/SiPixelTemplateDBObject.h" +#include "CondFormats/SiPixelTransient/interface/SiPixelTemplate.h" + #include #include @@ -52,9 +56,15 @@ class TrackingRecHitAlgorithm //this function will only be called once per stream virtual void beginStream(const edm::StreamID& id); + //this function will only be called once per run + virtual void beginRun(edm::Run const& run, const edm::EventSetup& eventSetup, + const SiPixelTemplateDBObject * pixelTemplateDBObjectPtr, + std::vector< SiPixelTemplateStore > & tempStoreRef ); + //this function will only be called once per event virtual void beginEvent(edm::Event& event, const edm::EventSetup& eventSetup); + //the main action is here virtual TrackingRecHitProductPtr process(TrackingRecHitProductPtr product) const; //this function will only be called once per event diff --git a/FastSimulation/TrackingRecHitProducer/plugins/PixelBarrelTemplateSmearerPlugin.cc b/FastSimulation/TrackingRecHitProducer/plugins/PixelBarrelTemplateSmearerPlugin.cc deleted file mode 100644 index 1bd32dd8f7e1b..0000000000000 --- a/FastSimulation/TrackingRecHitProducer/plugins/PixelBarrelTemplateSmearerPlugin.cc +++ /dev/null @@ -1,130 +0,0 @@ - -// SiPixel Gaussian Smearing -#include "FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h" -#include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithmFactory.h" -#include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitProduct.h" - -// Geometry -//#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" -#include "Geometry/CommonDetUnit/interface/GeomDet.h" -//#include "Geometry/CommonDetUnit/interface/GeomDetType.h" -#include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h" -#include "DataFormats/GeometryVector/interface/LocalPoint.h" -#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h" - -// Famos -#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h" -#include "FastSimulation/Utilities/interface/SimpleHistogramGenerator.h" - -// STL - -// ROOT -#include -#include -#include - - -using namespace std; - -class PixelBarrelTemplateSmearerPlugin: - public PixelTemplateSmearerBase -{ - public: - explicit PixelBarrelTemplateSmearerPlugin( - const std::string& name, - const edm::ParameterSet& config, - edm::ConsumesCollector& consumesCollector - ); - ~PixelBarrelTemplateSmearerPlugin() override; - - private: - void initializeBarrel(); -}; - - -PixelBarrelTemplateSmearerPlugin::PixelBarrelTemplateSmearerPlugin( - const std::string& name, - const edm::ParameterSet& config, - edm::ConsumesCollector& consumesCollector -): - PixelTemplateSmearerBase(name,config,consumesCollector) -{ - isForward = false; - - theBigPixelResolutionFileName = config.getParameter( "BigPixelBarrelResolutionFile" ); - theBigPixelResolutionFile = std::make_unique( edm::FileInPath( theBigPixelResolutionFileName ).fullPath().c_str() ,"READ"); - theEdgePixelResolutionFileName = config.getParameter( "EdgePixelBarrelResolutionFile" ); - theEdgePixelResolutionFile = std::make_unique( edm::FileInPath( theEdgePixelResolutionFileName ).fullPath().c_str() ,"READ"); - theRegularPixelResolutionFileName = config.getParameter( "RegularPixelBarrelResolutionFile" ); - theRegularPixelResolutionFile = std::make_unique( edm::FileInPath( theRegularPixelResolutionFileName ).fullPath().c_str() ,"READ"); - - theMergingProbabilityFileName = config.getParameter( "MergingProbabilityBarrelFile" ); - theMergingProbabilityFile = std::make_unique( edm::FileInPath( theMergingProbabilityFileName ).fullPath().c_str() ,"READ"); - theMergedPixelResolutionXFileName = config.getParameter( "MergedPixelBarrelResolutionXFile" ); - theMergedPixelResolutionXFile = std::make_unique( edm::FileInPath( theMergedPixelResolutionXFileName ).fullPath().c_str() ,"READ"); - theMergedPixelResolutionYFileName = config.getParameter( "MergedPixelBarrelResolutionYFile" ); - theMergedPixelResolutionYFile = std::make_unique( edm::FileInPath( theMergedPixelResolutionYFileName ).fullPath().c_str() ,"READ"); - - initializeBarrel(); - - if (!SiPixelTemplate::pushfile(templateId, thePixelTemp_)) - { - throw cms::Exception("PixelBarrelTemplateSmearerPlugin:")<<"SiPixel Barrel Template Not Loaded Correctly!"<Get(Form("DQMData/clustBPIX/hx%u",singleBigPixelHistN))); - theYHistos[singleBigPixelHistN] = new SimpleHistogramGenerator((TH1F*)theBigPixelResolutionFile->Get(Form("DQMData/clustBPIX/hy%u",singleBigPixelHistN))); - - unsigned int singlePixelHistN = 1*10000 + cotbetaHistBin*10 + cotalphaHistBin; - theXHistos[singlePixelHistN] = new SimpleHistogramGenerator((TH1F*)theRegularPixelResolutionFile->Get(Form("hx%u",singlePixelHistN))); - theYHistos[singlePixelHistN] = new SimpleHistogramGenerator((TH1F*)theRegularPixelResolutionFile->Get(Form("hy%u",singlePixelHistN))); - - for(unsigned qbinBin=1; qbinBin<=resqbin_binN; ++qbinBin ) - { - unsigned int edgePixelHistN = cotalphaHistBin*1000 + cotbetaHistBin*10 + qbinBin; - theXHistos[edgePixelHistN] = new SimpleHistogramGenerator((TH1F*)theEdgePixelResolutionFile->Get(Form("DQMData/clustBPIX/hx0%u",edgePixelHistN))); - theYHistos[edgePixelHistN] = new SimpleHistogramGenerator((TH1F*)theEdgePixelResolutionFile->Get(Form("DQMData/clustBPIX/hy0%u",edgePixelHistN))); - - unsigned int multiPixelBigHistN = 1*1000000 + 1*100000 + cotalphaHistBin*1000 + cotbetaHistBin * 10 + qbinBin; - theXHistos[multiPixelBigHistN] = new SimpleHistogramGenerator((TH1F*)theBigPixelResolutionFile->Get(Form("DQMData/clustBPIX/hx%u",multiPixelBigHistN))); - theYHistos[multiPixelBigHistN] = new SimpleHistogramGenerator((TH1F*)theBigPixelResolutionFile->Get(Form("DQMData/clustBPIX/hy%u",multiPixelBigHistN))); - - unsigned int multiPixelHistN = 1*100000 + 1*10000 + cotbetaHistBin*100 + cotalphaHistBin*10 + qbinBin; - theXHistos[multiPixelHistN] = new SimpleHistogramGenerator((TH1F*)theRegularPixelResolutionFile->Get(Form("hx%u",multiPixelHistN))); - theYHistos[multiPixelHistN] = new SimpleHistogramGenerator((TH1F*)theRegularPixelResolutionFile->Get(Form("hy%u",multiPixelHistN))); - } - } - } -} - - -DEFINE_EDM_PLUGIN( - TrackingRecHitAlgorithmFactory, - PixelBarrelTemplateSmearerPlugin, - "PixelBarrelTemplateSmearerPlugin" -); diff --git a/FastSimulation/TrackingRecHitProducer/plugins/PixelForwardTemplateSmearerPlugin.cc b/FastSimulation/TrackingRecHitProducer/plugins/PixelForwardTemplateSmearerPlugin.cc deleted file mode 100644 index 303f71db3e70a..0000000000000 --- a/FastSimulation/TrackingRecHitProducer/plugins/PixelForwardTemplateSmearerPlugin.cc +++ /dev/null @@ -1,150 +0,0 @@ - -/** PixelForwardTemplateSmearerPlugin.cc - * --------------------------------------------------------------------- - * Description: see PixelForwardTemplateSmearerPlugin.h - * Authors: R. Ranieri (CERN), M. Galanti - * History: Oct 11, 2006 - initial version - * - * New Pixel Resolution Parameterization - * Introduce SiPixelTemplate Object to Assign Pixel Errors - * by G. Hu - * --------------------------------------------------------------------- - */ - -// SiPixel Gaussian Smearing -#include "FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h" -#include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithmFactory.h" -#include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitProduct.h" - -// Geometry -//#include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" -#include "Geometry/CommonDetUnit/interface/GeomDet.h" -//#include "Geometry/CommonDetUnit/interface/GeomDetType.h" -#include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h" -#include "DataFormats/GeometryVector/interface/LocalPoint.h" -#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h" - -// Famos -#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h" -#include "FastSimulation/Utilities/interface/SimpleHistogramGenerator.h" - -// STL - -// ROOT -#include -#include -#include -//#include - -//#define FAMOS_DEBUG - -const double microntocm = 0.0001; -using namespace std; - - -class PixelForwardTemplateSmearerPlugin: - public PixelTemplateSmearerBase -{ - public: - explicit PixelForwardTemplateSmearerPlugin( - const std::string& name, - const edm::ParameterSet& config, - edm::ConsumesCollector& consumesCollector - ); - ~PixelForwardTemplateSmearerPlugin() override; - - private: - void initializeForward(); -}; - - -PixelForwardTemplateSmearerPlugin::PixelForwardTemplateSmearerPlugin( - const std::string& name, - const edm::ParameterSet& config, - edm::ConsumesCollector& consumesCollector -): - PixelTemplateSmearerBase(name,config,consumesCollector) -{ - - isForward = true; - - theEdgePixelResolutionFileName = config.getParameter( "EdgePixelForwardResolutionFile" ); - theEdgePixelResolutionFile = std::make_unique( edm::FileInPath( theEdgePixelResolutionFileName ).fullPath().c_str() ,"READ"); - theBigPixelResolutionFileName = config.getParameter( "BigPixelForwardResolutionFile" ); - theBigPixelResolutionFile = std::make_unique( edm::FileInPath( theBigPixelResolutionFileName ).fullPath().c_str() ,"READ"); - theRegularPixelResolutionFileName = config.getParameter( "RegularPixelForwardResolutionFile" ); - theRegularPixelResolutionFile = std::make_unique( edm::FileInPath( theRegularPixelResolutionFileName ).fullPath().c_str() ,"READ"); - - theMergingProbabilityFileName = config.getParameter( "MergingProbabilityForwardFile" ); - theMergingProbabilityFile =std::make_unique( edm::FileInPath( theMergingProbabilityFileName ).fullPath().c_str() ,"READ"); - theMergedPixelResolutionXFileName = config.getParameter( "MergedPixelForwardResolutionXFile" ); - theMergedPixelResolutionXFile = std::make_unique( edm::FileInPath( theMergedPixelResolutionXFileName ).fullPath().c_str() ,"READ"); - theMergedPixelResolutionYFileName = config.getParameter( "MergedPixelForwardResolutionYFile" ); - theMergedPixelResolutionYFile = std::make_unique( edm::FileInPath( theMergedPixelResolutionYFileName ).fullPath().c_str() ,"READ"); - - initializeForward(); - - - - if( ! SiPixelTemplate::pushfile(templateId, thePixelTemp_) ) - { - throw cms::Exception("PixelForwardTemplateSmearerPlugin:") <<"SiPixel Forward Template Not Loaded Correctly!"; - } -} - - -PixelForwardTemplateSmearerPlugin::~PixelForwardTemplateSmearerPlugin() -{ -} - - -void PixelForwardTemplateSmearerPlugin::initializeForward() -{ - rescotAlpha_binMin = 0.1; - rescotAlpha_binWidth = 0.1; - rescotAlpha_binN = 4; - rescotBeta_binMin = 0.; - rescotBeta_binWidth = 0.15; - rescotBeta_binN = 4; - resqbin_binMin = 0; - resqbin_binWidth = 1; - resqbin_binN = 4; - - // Initialize the forward histos once and for all, and prepare the random generation - for ( unsigned cotalphaHistBin=1; cotalphaHistBin<=rescotAlpha_binN; ++cotalphaHistBin) - { - for ( unsigned cotbetaHistBin=1; cotbetaHistBin<=rescotBeta_binN; ++cotbetaHistBin) - { - for( unsigned qbinBin=1; qbinBin<=resqbin_binN; ++qbinBin ) - { - unsigned int edgePixelHistN = cotalphaHistBin*1000 + cotbetaHistBin*10 + qbinBin; - theXHistos[edgePixelHistN] = new SimpleHistogramGenerator((TH1F*)theEdgePixelResolutionFile->Get(Form("DQMData/clustFPIX/fhx0%u",edgePixelHistN))); - theYHistos[edgePixelHistN] = new SimpleHistogramGenerator((TH1F*)theEdgePixelResolutionFile->Get(Form("DQMData/clustFPIX/fhy0%u",edgePixelHistN))); - - unsigned int PixelHistN = 10000 + cotbetaHistBin*100 + cotalphaHistBin*10 + qbinBin; - theXHistos[PixelHistN] = new SimpleHistogramGenerator((TH1F*) theRegularPixelResolutionFile->Get(Form("hx0%u",PixelHistN))); - theYHistos[PixelHistN] = new SimpleHistogramGenerator((TH1F*) theRegularPixelResolutionFile->Get(Form("hy0%u",PixelHistN))); - } - } - } - - for ( unsigned cotalphaHistBin=1; cotalphaHistBin<=rescotAlpha_binN; ++cotalphaHistBin) - { - for ( unsigned cotbetaHistBin=1; cotbetaHistBin<=rescotBeta_binN; ++cotbetaHistBin) - { - unsigned int SingleBigPixelHistN = 100000 + cotalphaHistBin*100 + cotbetaHistBin; - theXHistos[SingleBigPixelHistN] = new SimpleHistogramGenerator((TH1F*)theBigPixelResolutionFile->Get(Form("DQMData/clustFPIX/fhx%u",SingleBigPixelHistN))); - theYHistos[SingleBigPixelHistN] = new SimpleHistogramGenerator((TH1F*)theBigPixelResolutionFile->Get(Form("DQMData/clustFPIX/fhy%u",SingleBigPixelHistN))); - - unsigned int SinglePixelHistN = cotbetaHistBin*10 + cotalphaHistBin; - theXHistos[SinglePixelHistN] = new SimpleHistogramGenerator((TH1F*)theRegularPixelResolutionFile->Get(Form("hx000%u",SinglePixelHistN))); - theYHistos[SinglePixelHistN] = new SimpleHistogramGenerator((TH1F*)theRegularPixelResolutionFile->Get(Form("hy000%u",SinglePixelHistN))); - } - } -} - -DEFINE_EDM_PLUGIN( - TrackingRecHitAlgorithmFactory, - PixelForwardTemplateSmearerPlugin, - "PixelForwardTemplateSmearerPlugin" -); diff --git a/FastSimulation/TrackingRecHitProducer/plugins/PixelTemplateSmearerPlugin.cc b/FastSimulation/TrackingRecHitProducer/plugins/PixelTemplateSmearerPlugin.cc new file mode 100644 index 0000000000000..2ac1151120763 --- /dev/null +++ b/FastSimulation/TrackingRecHitProducer/plugins/PixelTemplateSmearerPlugin.cc @@ -0,0 +1,51 @@ + +// SiPixel Gaussian Smearing +#include "FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h" +#include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithmFactory.h" +#include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitProduct.h" +#include "FastSimulation/TrackingRecHitProducer/interface/PixelResolutionHistograms.h" + +// Geometry +#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h" +/// If we ever need to port back to 9X: #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h" +#include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h" +#include "DataFormats/GeometryVector/interface/LocalPoint.h" +#include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h" + +// Famos +#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h" +#include "FastSimulation/Utilities/interface/SimpleHistogramGenerator.h" + + +class PixelTemplateSmearerPlugin: + public PixelTemplateSmearerBase +{ +public: + explicit PixelTemplateSmearerPlugin( const std::string& name, + const edm::ParameterSet& config, + edm::ConsumesCollector& consumesCollector + ); + ~PixelTemplateSmearerPlugin() override; +}; + + +PixelTemplateSmearerPlugin::PixelTemplateSmearerPlugin( + const std::string& name, + const edm::ParameterSet& config, + edm::ConsumesCollector& consumesCollector +): + PixelTemplateSmearerBase(name, config, consumesCollector) +{ +} + + +PixelTemplateSmearerPlugin::~PixelTemplateSmearerPlugin() +{ +} + + +DEFINE_EDM_PLUGIN( + TrackingRecHitAlgorithmFactory, + PixelTemplateSmearerPlugin, + "PixelTemplateSmearerPlugin" +); diff --git a/FastSimulation/TrackingRecHitProducer/plugins/TrackingRecHitProducer.cc b/FastSimulation/TrackingRecHitProducer/plugins/TrackingRecHitProducer.cc index 940d6835c252e..be9c780161648 100644 --- a/FastSimulation/TrackingRecHitProducer/plugins/TrackingRecHitProducer.cc +++ b/FastSimulation/TrackingRecHitProducer/plugins/TrackingRecHitProducer.cc @@ -13,6 +13,7 @@ #include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithm.h" #include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithmFactory.h" #include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitPipe.h" +#include "FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h" #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h" #include "DataFormats/TrackerCommon/interface/TrackerTopology.h" @@ -27,6 +28,10 @@ #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h" #include "Geometry/Records/interface/TrackerTopologyRcd.h" +// Pixel-related stuff: +#include "CondFormats/SiPixelObjects/interface/SiPixelTemplateDBObject.h" +#include "CondFormats/SiPixelTransient/interface/SiPixelTemplate.h" +#include "CalibTracker/Records/interface/SiPixelTemplateDBObjectESProducerRcd.h" #include "FastSimulation/TrackingRecHitProducer/interface/TrackerDetIdSelector.h" @@ -46,6 +51,7 @@ class TrackingRecHitProducer: edm::IOVSyncValue _iovSyncValue; std::map _detIdPipes; void setupDetIdPipes(const edm::EventSetup& eventSetup); + std::vector< SiPixelTemplateStore > _pixelTempStore ; // pixel template storage public: TrackingRecHitProducer(const edm::ParameterSet& config); @@ -99,8 +105,12 @@ TrackingRecHitProducer::~TrackingRecHitProducer() delete algo; } _recHitAlgorithms.clear(); + + //--- Delete the templates. This is safe even if thePixelTemp_ vector is empty. + for (auto x : _pixelTempStore) x.destroy(); } + void TrackingRecHitProducer::beginStream(edm::StreamID id) { for (TrackingRecHitAlgorithm* algo: _recHitAlgorithms) @@ -109,9 +119,30 @@ void TrackingRecHitProducer::beginStream(edm::StreamID id) } } -void TrackingRecHitProducer::beginRun(edm::Run const&, const edm::EventSetup& eventSetup) +void TrackingRecHitProducer::beginRun(edm::Run const& run, const edm::EventSetup& eventSetup) { + //--- Since all pixel algorithms (of which there could be several) use the same + // templateStore, filled out from the same DB Object, we need to it centrally + // (namely here), and then distribute it to the algorithms. Note that only + // the pixel algorithms implement beginRun(), for the strip tracker this defaults + // to a no-op. + + edm::ESHandle templateDBobject; + eventSetup.get().get(templateDBobject); + const SiPixelTemplateDBObject * pixelTemplateDBObject = templateDBobject.product(); + + //--- Now that we have the DB object, load the correct templates from the DB. + // (They are needed for data and full sim MC, so in a production FastSim + // run, everything should already be in the DB.) + if ( !SiPixelTemplate::pushfile( *pixelTemplateDBObject, _pixelTempStore ) ) { + throw cms::Exception("TrackingRecHitProducer:") + << "SiPixel Templates not loaded correctly from the DB object!" << std::endl; + } + for (TrackingRecHitAlgorithm* algo: _recHitAlgorithms) + { + algo->beginRun(run, eventSetup, pixelTemplateDBObject, _pixelTempStore ); + } } void TrackingRecHitProducer::setupDetIdPipes(const edm::EventSetup& eventSetup) diff --git a/FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase0_cfi.py b/FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase0_cfi.py new file mode 100644 index 0000000000000..48e71e84c70df --- /dev/null +++ b/FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase0_cfi.py @@ -0,0 +1,59 @@ +import FWCore.ParameterSet.Config as cms + + +pixelPluginsPhase0=cms.VPSet() + +#===================================================================================== +#--- Phase 0 Pixel Barrel +# +# Layer Template Cluster file Resolution histograms +# ----------------------------------------------------------------------------- +# * 5611 template_events_d39921.out.gz pixel_histos39921_5611_6.root +# +# +pixelPluginsPhase0.append( + cms.PSet( + select = cms.string("subdetId==BPX"), + isBarrel = cms.bool(True), + name = cms.string("pixelSmearerBarrel"), + type = cms.string("PixelTemplateSmearerPlugin"), + # templateId = cms.int32( 5611 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos39921_5611_6.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelBig2017.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelEdge2017.root'), + # + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bysmear.root'), + ) +) + + +#===================================================================================== +#--- Phase 0 Pixel Forward +# +# Panel Template Cluster file Resolution histograms +# ----------------------------------------------------------------------------- +# * 6415 template_events_d40722.out.gz pixel_histos40722_6415_6.root +# +# +pixelPluginsPhase0.append( + cms.PSet( + select=cms.string("subdetId==FPX"), + isBarrel = cms.bool(False), + name = cms.string("pixelSmearerForward"), + type = cms.string("PixelTemplateSmearerPlugin"), + # templateId = cms.int32( 6415 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos40722_6415_6.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/ForwardBig2017.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/ForwardEdge2017.root'), + # + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fysmear.root'), + ) +) + + diff --git a/FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase1_cfi.py b/FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase1_cfi.py new file mode 100644 index 0000000000000..d000b161f8a34 --- /dev/null +++ b/FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase1_cfi.py @@ -0,0 +1,191 @@ +import FWCore.ParameterSet.Config as cms + +pixelPluginsPhase1=cms.VPSet() + + +#===================================================================================== +#--- Phase 1 Pixel Barrel +# +# Layer Template Cluster file Resolution histograms +# ----------------------------------------------------------------------------- +# BPL1 2403 template_events_d63207.out.gz pixel_histos63207_2403.root +# +#--- Layer 1 +# +pixelPluginsPhase1.append( + cms.PSet( + select = cms.string("subdetId==BPX && pxbLayer==1"), + isBarrel = cms.bool(True), + name = cms.string("pixelSmearerBarrelLayer1"), + type = cms.string("PixelTemplateSmearerPlugin"), + # templateId = cms.int32( 2403 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos63207_2403_6.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelBig2017.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelEdge2017.root'), + # + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bysmear.root'), + ) +) + +# +#--- Layer 2 +# BPL2 2413 template_events_d63507.out.gz pixel_histos63507_2413.root +# +pixelPluginsPhase1.append( + cms.PSet( + select = cms.string("subdetId==BPX && pxbLayer==2"), + isBarrel = cms.bool(True), + name = cms.string("pixelSmearerBarrelLayer2"), + type = cms.string("PixelTemplateSmearerPlugin"), + # templateId = cms.int32( 2413 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos63507_2413_6.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelBig2017.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelEdge2017.root'), + # + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bysmear.root'), + ) +) + +# +#--- Layer 3 +# BPL3 2423 template_events_d63807.out.gz pixel_histos63807_2423.root +# +pixelPluginsPhase1.append( + cms.PSet( + select = cms.string("subdetId==BPX && pxbLayer==3"), + isBarrel = cms.bool(True), + name = cms.string("pixelSmearerBarrelLayer3"), + type = cms.string("PixelTemplateSmearerPlugin"), + # templateId = cms.int32( 2413 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos63807_2423_6.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelBig2017.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelEdge2017.root'), + # + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bysmear.root'), + ) +) + + +# +#--- Layer 4 +# BPL4 2433 template_events_d63807.out.gz pixel_histos64107_2433.root +# +pixelPluginsPhase1.append( + cms.PSet( + select = cms.string("subdetId==BPX && pxbLayer==4"), + isBarrel = cms.bool(True), + name = cms.string("pixelSmearerBarrelLayer4"), + type = cms.string("PixelTemplateSmearerPlugin"), + # templateId = cms.int32( 2413 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos64107_2433_6.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelBig2017.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelEdge2017.root'), + # + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bysmear.root'), + ) +) + + + + +#===================================================================================== +#--- Phase 1 Pixel Forward +# +# Panel Template Cluster file Resolution histograms +# ----------------------------------------------------------------------------- +# FPR2P1 2443 template_events_d64237.out.gz pixel_histos64237_2443.root +# +#--- Ring 2, Panel 1 +pixelPluginsPhase1.append( + cms.PSet( + select=cms.string("subdetId==FPX && pxfBlade>22 && pxfPanel==1"), ## 1-56 (Ring 1 is 1-22, Ring 2 is 23-56) + isBarrel = cms.bool(False), + name = cms.string("pixelSmearerForwardRing2Panel1"), + type = cms.string("PixelTemplateSmearerPlugin"), + # templateId = cms.int32( 2443 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos64237_2443_6.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/ForwardBig2017.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/ForwardEdge2017.root'), + # + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fysmear.root'), + ) +) + + + +#--- Ring 1, Panel 1 +# FPR1P1 2453 template_events_d64367.out.gz pixel_histos64367_2453.root +pixelPluginsPhase1.append( + cms.PSet( + select=cms.string("subdetId==FPX && pxfBlade<23 && pxfPanel==1"), ## 1-56 (Ring 1 is 1-22, Ring 2 is 23-56) + isBarrel = cms.bool(False), + name = cms.string("pixelSmearerForwardRing1Panel1"), + type = cms.string("PixelTemplateSmearerPlugin"), + # templateId = cms.int32( 2453 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos64367_2453_6.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/ForwardBig2017.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/ForwardEdge2017.root'), + # + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fysmear.root'), + ) +) + + +#--- Ring 1, Panel 2 +# FPR1P2 2463 template_events_d64497.out.gz pixel_histos64497_2463.root +pixelPluginsPhase1.append( + cms.PSet( + select=cms.string("subdetId==FPX && pxfBlade<23 && pxfPanel==2"), ## 1-56 (Ring 1 is 1-22, Ring 2 is 23-56) + isBarrel = cms.bool(False), + name = cms.string("pixelSmearerForwardRing1Panel2"), + type = cms.string("PixelTemplateSmearerPlugin"), + # templateId = cms.int32( 2463 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos64497_2463_6.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/ForwardBig2017.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/ForwardEdge2017.root'), + # + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fysmear.root'), + ) +) + +#--- Ring 2, Panel 2 +# FPR2P2 2473 template_events_d64627.out.gz pixel_histos64627_2473.root +pixelPluginsPhase1.append( + cms.PSet( + select=cms.string("subdetId==FPX && pxfBlade>22 && pxfPanel==2"), ## 1-56 (Ring 1 is 1-22, Ring 2 is 23-56) + isBarrel = cms.bool(False), + name = cms.string("pixelSmearerForwardRing2Panel2"), + type = cms.string("PixelTemplateSmearerPlugin"), + # templateId = cms.int32( 2473 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos64627_2473_6.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/ForwardBig2017.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/ForwardEdge2017.root'), + # + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fysmear.root'), + ) +) + diff --git a/FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase2_cfi.py b/FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase2_cfi.py new file mode 100644 index 0000000000000..81c23ff62073a --- /dev/null +++ b/FastSimulation/TrackingRecHitProducer/python/PixelPluginsPhase2_cfi.py @@ -0,0 +1,46 @@ +import FWCore.ParameterSet.Config as cms + + +pixelPluginsPhase2=cms.VPSet() + + +pixelPluginsPhase2.append( + cms.PSet( + select = cms.string("subdetId==BPX"), + isBarrel = cms.bool(True), + name = cms.string("pixelBarrelSmearer"), + type = cms.string("PixelTemplateSmearerPlugin"), + templateId = cms.int32( 40 ), + ## templateId = cms.int32( 292 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos58360_292.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos58360_292.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos58360_292.root'), + ## BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelEdge2017.root'), + ## EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelEdge2017.root'), + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bysmear.root'), + ) +) + +pixelPluginsPhase2.append( + cms.PSet( + select=cms.string("subdetId==FPX"), + isBarrel = cms.bool(False), + name = cms.string("pixelForwardSmearer"), + type = cms.string("PixelTemplateSmearerPlugin"), + templateId = cms.int32( 41 ), + ## templateId = cms.int32( 291 ), + RegularPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos49765_291.root'), + BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos49765_291.root'), + EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/pixel_histos49765_291.root'), + ## BigPixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelEdge2017.root'), + ## EdgePixelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/BarrelEdge2017.root'), + MergeHitsOn = cms.bool(False), + MergingProbabilityFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fmergeprob.root'), + MergedPixelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fxsmear.root'), + MergedPixelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fysmear.root'), + ) +) + diff --git a/FastSimulation/TrackingRecHitProducer/python/TrackingRecHitProducer_cfi.py b/FastSimulation/TrackingRecHitProducer/python/TrackingRecHitProducer_cfi.py index 19a3996896b83..ddf4ede192576 100644 --- a/FastSimulation/TrackingRecHitProducer/python/TrackingRecHitProducer_cfi.py +++ b/FastSimulation/TrackingRecHitProducer/python/TrackingRecHitProducer_cfi.py @@ -1,57 +1,33 @@ +# Python 2 vs 3 compatibility library: +import six + import FWCore.ParameterSet.Config as cms +# Load the detailed configurations of pixel plugins. +# NB: for any new detector geometry (e.g. Phase 2 varians), we should write a new plugin +# config file, and import it here, and below use its own Era to load it. +# +from FastSimulation.TrackingRecHitProducer.PixelPluginsPhase0_cfi import pixelPluginsPhase0 +from FastSimulation.TrackingRecHitProducer.PixelPluginsPhase1_cfi import pixelPluginsPhase1 +from FastSimulation.TrackingRecHitProducer.PixelPluginsPhase2_cfi import pixelPluginsPhase2 + +# The default is (for better of worse) Phase 0: +# fastTrackerRecHits = cms.EDProducer("TrackingRecHitProducer", simHits = cms.InputTag("fastSimProducer","TrackerHits"), - plugins=cms.VPSet() + plugins = pixelPluginsPhase0 ) -fastTrackerRecHits.plugins.append( - cms.PSet( - name = cms.string("pixelBarrelSmearer"), - type=cms.string("PixelBarrelTemplateSmearerPlugin"), - BigPixelBarrelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/NewPixelResolutionBarrel38T.root'), - EdgePixelBarrelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/NewPixelResolutionBarrelEdge38T.root'), - RegularPixelBarrelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/PixelBarrelResolution2014.root'), - BigPixelForwardResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/NewPixelResolutionForward38T.root'), - EdgePixelForwardResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/NewPixelResolutionForward38T.root'), - RegularPixelForwardResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/PixelForwardResolution2014.root'), - UseCMSSWPixelParametrization = cms.bool(False), - MergeHitsOn = cms.bool(False), - MergingProbabilityBarrelFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bmergeprob.root'), - MergingProbabilityForwardFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fmergeprob.root'), - MergedPixelBarrelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bxsmear.root'), - MergedPixelForwardResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fxsmear.root'), - MergedPixelBarrelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bysmear.root'), - MergedPixelForwardResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fysmear.root'), - templateId = cms.int32( 40 ), - select=cms.string("subdetId==BPX"), - ) -) +# Phase 1 Era: replace plugins by Phase 1 plugins +from Configuration.Eras.Modifier_phase1Pixel_cff import phase1Pixel +phase1Pixel.toModify(fastTrackerRecHits, plugins = pixelPluginsPhase1) -fastTrackerRecHits.plugins.append( - cms.PSet( - name = cms.string("pixelForwardSmearer"), - type=cms.string("PixelForwardTemplateSmearerPlugin"), - BigPixelBarrelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/NewPixelResolutionBarrel38T.root'), - EdgePixelBarrelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/NewPixelResolutionBarrelEdge38T.root'), - RegularPixelBarrelResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/PixelBarrelResolution2014.root'), - BigPixelForwardResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/NewPixelResolutionForward38T.root'), - EdgePixelForwardResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/NewPixelResolutionForward38T.root'), - RegularPixelForwardResolutionFile = cms.string('FastSimulation/TrackingRecHitProducer/data/PixelForwardResolution2014.root'), - UseCMSSWPixelParametrization = cms.bool(False), - MergeHitsOn = cms.bool(False), - MergingProbabilityBarrelFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bmergeprob.root'), - MergingProbabilityForwardFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fmergeprob.root'), - MergedPixelBarrelResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bxsmear.root'), - MergedPixelForwardResolutionXFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fxsmear.root'), - MergedPixelBarrelResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/bysmear.root'), - MergedPixelForwardResolutionYFile = cms.string('FastSimulation/TrackingRecHitProducer/data/fysmear.root'), - templateId = cms.int32( 41 ), - select=cms.string("subdetId==FPX"), - ) -) +# Phase 2 Era: replace plugins by Phase 2 plugins, etc... +from Configuration.Eras.Modifier_phase2_tracker_cff import phase2_tracker +phase2_tracker.toModify(fastTrackerRecHits, plugins = pixelPluginsPhase2) +# Configure strip tracker Gaussian-smearing plugins: trackerStripGaussianResolutions={ "TIB": { 1: cms.double(0.00195), diff --git a/FastSimulation/TrackingRecHitProducer/src/PixelResolutionHistograms.cc b/FastSimulation/TrackingRecHitProducer/src/PixelResolutionHistograms.cc new file mode 100644 index 0000000000000..0424a80df0a7a --- /dev/null +++ b/FastSimulation/TrackingRecHitProducer/src/PixelResolutionHistograms.cc @@ -0,0 +1,580 @@ +// +// The implementation of the PixelResolutinHistograms.cc class. Please +// look at PixelResolutionHistograms.h header file for the interface. +// +//------------------------------------------------------------------------------ + +// The switch, undefined in CMSSW release, and defined by standalone compilation script: + + +#ifdef SI_PIXEL_TEMPLATE_STANDALONE +// +//--- Stand-alone: Include a the header file from the local directory, as well as +// dummy implementations of SimpleHistogramGenerator, LogInfo, LogError and LogDebug... +// +#include "PixelResolutionHistograms.h" +// +class TH1F; +class TH2F; +class SimpleHistogramGenerator { +public: + SimpleHistogramGenerator(TH1F * hist) : hist_(hist) {}; +private: + TH1F * hist_; // we don't own it +}; +#define LOGDEBUG std::cout +#define LOGERROR std::cout +#define LOGINFO std::cout +// +#else +//--- We're inside a CMSSW release: Include the real thing. +// +#include "FastSimulation/TrackingRecHitProducer/interface/PixelResolutionHistograms.h" +#include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h" +#include "FastSimulation/Utilities/interface/SimpleHistogramGenerator.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +// +#define LOGDEBUG LogDebug("") +#define LOGERROR edm::LogError("Error") +#define LOGINFO edm::LogInfo("Info") +// +#endif + + + + +// Generic C stuff +#include +#include +#include + +// ROOT +#include +#include +#include + +// Global definitions +const float cmtomicron = 10000.0; + + +//------------------------------------------------------------------------------ +// Constructor: Books the FastSim Histograms, given the input parameters +// which are provided as arguments. These variables are then const inside +// the class. (That is, once we make the histograms, we can't change the +// definition of the binning.) +//------------------------------------------------------------------------------ +PixelResolutionHistograms:: +PixelResolutionHistograms( std::string filename, // ROOT file for histograms + std::string rootdir, // Subdirectory in the file, "" if none + std::string descTitle, // Descriptive title + unsigned int detType, // Where we are... (&&& do we need this?) + double cotbetaBinWidth, + double cotbetaLowEdge, + int cotbetaBins, + double cotalphaBinWidth, + double cotalphaLowEdge, + int cotalphaBins) //, + //int qbinWidth, + //int qbins ) + : weOwnHistograms_(true), // we'll be making some histos + detType_ ( detType ), + cotbetaBinWidth_ ( cotbetaBinWidth ), + cotbetaLowEdge_ ( cotbetaLowEdge ), + cotbetaBins_ ( cotbetaBins ), + cotalphaBinWidth_ ( cotalphaBinWidth ), + cotalphaLowEdge_ ( cotalphaLowEdge ), + cotalphaBins_ ( cotalphaBins ), + qbinWidth_ ( 1 ), + qbins_ ( 4 ), + binningHisto_(nullptr), + resMultiPixelXHist_(), resSinglePixelXHist_(), // all to nullptr + resMultiPixelYHist_(), resSinglePixelYHist_(), // all to nullptr + qbinHist_(), // all to nullptr + file_(nullptr), + status_(0), + resMultiPixelXGen_(), resSinglePixelXGen_(), + resMultiPixelYGen_(), resSinglePixelYGen_(), + qbinGen_() +{ + + file_ = std::make_unique( filename.c_str(), "RECREATE" ); + //Resolution binning + // const double cotbetaBinWidth = 1.0; + // const double cotbetaLowEdge = -11.5 ; + // const int cotbetaBins = 23; + // const double cotalphaBinWidth = 0.08 ; + // const double cotalphaLowEdge = -0.36 ; + // const int cotalphaBins = 9; + // const int qbinWidth = 1; + // const int qbins = 4; + + // Dummy 2D histogram to store binning: + binningHisto_ = new TH2F( "ResHistoBinning", descTitle.c_str(), + cotbetaBins, cotbetaLowEdge, cotbetaLowEdge + cotbetaBins*cotbetaBinWidth, + cotalphaBins, cotalphaLowEdge, cotalphaLowEdge + cotalphaBins*cotalphaBinWidth ); + + // Store detType in the underflow bin + binningHisto_->SetBinContent(0, 0, detType_); + + // All other histograms: + Char_t histo[200]; + Char_t title[200]; + // + //--- Histograms for clusters with multiple pixels hit in a given direction. + // + for( int ii=0; ii1 X", + cotbetaLowEdge_ + ii*cotbetaBinWidth_ , cotbetaLowEdge_ + (ii+1)*cotbetaBinWidth_, + cotalphaLowEdge_ +jj*cotalphaBinWidth_, cotalphaLowEdge_ +(jj+1)*cotalphaBinWidth_, + kk+1 ); + // + resMultiPixelXHist_ [ ii ][ jj ][ kk ] = new TH1F(histo, title, 1000, -0.05, 0.05); + + sprintf( histo, "hy%d1%02d%d%d", detType_, ii+1, jj+1, kk+1 ); + sprintf( title, "cotbeta %.1f-%.1f cotalpha %.2f-%.2f qbin %d npixel>1 Y", + cotbetaLowEdge_ + ii*cotbetaBinWidth_ , cotbetaLowEdge_ + (ii+1)*cotbetaBinWidth_, + cotalphaLowEdge_ +jj*cotalphaBinWidth_, cotalphaLowEdge_ +(jj+1)*cotalphaBinWidth_, + kk+1 ); + // + resMultiPixelYHist_ [ ii ][ jj ][ kk ] = new TH1F(histo, title, 1000, -0.05, 0.05); + } + } + } + + + // + //--- Histograms for clusters where only a single pixel was hit in a given direction. + // + for( int ii=0; ii( filename.c_str() ,"READ"); + if ( !file_ ) { + status_ = 1; + LOGERROR << "PixelResolutionHistograms:: Error, file " << filename << " not found."; + return; // PixelTemplateSmearerBase will throw an exception upon our return. + } + + //--- The dummy 2D histogram with the binning of cot\beta and cot\alpha: + binningHisto_ = (TH2F*) file_->Get( Form( "%s%s" , rootdir.c_str(), "ResHistoBinning" ) ); + if ( !binningHisto_ ) { + status_ = 11; + LOGERROR << "PixelResolutionHistograms:: Error, binning histogrram ResHistoBinning not found."; + return; // PixelTemplateSmearerBase will throw an exception upon our return. + } + + if ( detType == -1 ) { + //--- Fish out detType from the underflow bin: + detType_ = binningHisto_->GetBinContent(0, 0); + } + else { + detType_ = detType; // constructor's argument overrides what's in ResHistoBinning histogram. + } + + //--- Now we fill the binning variables: + cotbetaAxis_ = binningHisto_->GetXaxis(); + cotbetaBinWidth_ = binningHisto_->GetXaxis()->GetBinWidth(1); // assume all same width + cotbetaLowEdge_ = binningHisto_->GetXaxis()->GetXmin(); // low edge of the first bin + cotbetaBins_ = binningHisto_->GetXaxis()->GetNbins(); + cotalphaAxis_ = binningHisto_->GetYaxis(); + cotalphaBinWidth_ = binningHisto_->GetYaxis()->GetBinWidth(1); // assume all same width; + cotalphaLowEdge_ = binningHisto_->GetYaxis()->GetXmin(); // low edge of the first bin; + cotalphaBins_ = binningHisto_->GetYaxis()->GetNbins(); + + //--- We want the following information to show up in *every* log file! + // LOGINFO << std::endl + std::cout << std::endl + << "Loading pixel resolution file = " << filename << std::endl + << " cotBeta[" << cotbetaLowEdge_ <<","<< cotbetaBinWidth_ <<","<< cotbetaBins_ << "]" << std::endl + << " cotAlpha[" << cotalphaLowEdge_ <<","<< cotalphaBinWidth_ <<","<< cotalphaBins_ << "]" + << std::endl; + + + if ( !ignore_multi ) { + // + //--- Histograms for clusters with multiple pixels hit in a given direction. + // + for( int ii=0; ii1 X", + cotbetaLowEdge_ + ii*cotbetaBinWidth_ , cotbetaLowEdge_ + (ii+1)*cotbetaBinWidth_, + cotalphaLowEdge_ +jj*cotalphaBinWidth_, cotalphaLowEdge_ +(jj+1)*cotalphaBinWidth_, + kk+1 ); + // + tmphist = (TH1F*) file_->Get( Form( "%s%s" , rootdir.c_str(), histo ) ); + if ( !tmphist ) { + status_ = 2; + LOGERROR << "Failed to find histogram=" << std::string( histo ); + return; + } + LOGDEBUG << "Found histo " << std::string(histo) + << " with title = " << std::string( tmphist->GetTitle() ) << std::endl; + if ( tmphist->GetEntries() < 5 ) { + LOGINFO << "Histogram " << std::string(histo) << " has only " << tmphist->GetEntries() + << " entries. Trouble ahead." << std::endl; + } + resMultiPixelXHist_ [ ii ][ jj ][ kk ] = tmphist; + resMultiPixelXGen_ [ ii ][ jj ][ kk ] = new SimpleHistogramGenerator( tmphist ); + + + sprintf( histo, "hy%d1%02d%d%d", detType_, ii+1, jj+1, kk+1 ); + sprintf( title, "cotbeta %.1f-%.1f cotalpha %.2f-%.2f qbin %d npixel>1 Y", + cotbetaLowEdge_ + ii*cotbetaBinWidth_ , cotbetaLowEdge_ + (ii+1)*cotbetaBinWidth_, + cotalphaLowEdge_ +jj*cotalphaBinWidth_, cotalphaLowEdge_ +(jj+1)*cotalphaBinWidth_, + kk+1 ); + // + tmphist = (TH1F*) file_->Get( Form( "%s%s" , rootdir.c_str(), histo ) ); + if ( !tmphist ) { + status_ = 3; + LOGERROR << "Failed to find histogram=" << std::string( histo ); + return; + } + LOGDEBUG << "Found histo " << std::string(histo) + << " with title = " << std::string( tmphist->GetTitle() ) << std::endl; + if ( tmphist->GetEntries() < 5 ) { + LOGINFO << "Histogram " << std::string(histo) << " has only " << tmphist->GetEntries() + << " entries. Trouble ahead." << std::endl; + } + resMultiPixelYHist_ [ ii ][ jj ][ kk ] = tmphist; + resMultiPixelYGen_ [ ii ][ jj ][ kk ] = new SimpleHistogramGenerator( tmphist ); + } + } + } + // + } // if (not ignore multi) + + + // + //--- Histograms for clusters where only a single pixel was hit in a given direction. + // + for( int ii=0; iiGet( Form( "%s%s" , rootdir.c_str(), histo ) ); + if ( !tmphist ) { + if ( !ignore_single ) { + LOGERROR << "Failed to find histogram=" << std::string( histo ); + status_ = 4; + return; + } + } + else { + LOGDEBUG << "Found histo " << std::string(histo) + << " with title = " << std::string( tmphist->GetTitle() ) << std::endl; + LOGDEBUG << "Found histo with title = " << std::string( tmphist->GetTitle() ) << std::endl; + if ( tmphist->GetEntries() < 5 ) { + LOGINFO << "Histogram " << std::string(histo) << " has only " << tmphist->GetEntries() + << " entries. Trouble ahead." << std::endl; + } + resSinglePixelXHist_ [ ii ][ jj ] = tmphist; + resSinglePixelXGen_ [ ii ][ jj ] = new SimpleHistogramGenerator( tmphist ); + } + + + //--- Single pixel, along Y. + // + sprintf( histo, "hy%d0%02d%d", detType_, ii+1, jj+1 ); + sprintf( title, "cotbeta %.1f-%.1f cotalpha %.2f-%.2f npixel=1 Y", + cotbetaLowEdge_ + ii*cotbetaBinWidth_ , cotbetaLowEdge_ + (ii+1)*cotbetaBinWidth_, + cotalphaLowEdge_ +jj*cotalphaBinWidth_, cotalphaLowEdge_ +(jj+1)*cotalphaBinWidth_ ); + // + tmphist = (TH1F*) file_->Get( Form( "%s%s" , rootdir.c_str(), histo ) ); + if ( !tmphist ) { + if ( !ignore_single ) { + LOGERROR << "Failed to find histogram=" << std::string( histo ); + status_ = 5; + return; + } + } + else { + LOGDEBUG << "Found histo " << std::string(histo) + << " with title = " << std::string( tmphist->GetTitle() ) << std::endl; + if ( tmphist->GetEntries() < 5 ) { + LOGINFO << "Histogram " << std::string(histo) << " has only " << tmphist->GetEntries() + << " entries. Trouble ahead." << std::endl; + } + resSinglePixelYHist_ [ ii ][ jj ] = tmphist; + resSinglePixelYGen_ [ ii ][ jj ] = new SimpleHistogramGenerator( tmphist ); + } + + + //--- qBin distribution, for this (cotbeta, cotalpha) bin. + // + sprintf( histo, "hqbin%d%02d%d", detType_, ii+1, jj+1 ); + sprintf( title, "cotbeta %.1f-%.1f cotalpha %.2f-%.2f qbin", + cotbetaLowEdge_ + ii*cotbetaBinWidth_ , cotbetaLowEdge_ + (ii+1)*cotbetaBinWidth_, + cotalphaLowEdge_ +jj*cotalphaBinWidth_, cotalphaLowEdge_ +(jj+1)*cotalphaBinWidth_ ); + // + tmphist = (TH1F*) file_->Get( Form( "%s%s" , rootdir.c_str(), histo ) ); + if ( !tmphist ) { + if ( !ignore_qBin ) { + LOGERROR << "Failed to find histogram=" << std::string( histo ); + status_ = 6; + return; + } + } + else { + LOGDEBUG << "Found histo " << std::string(histo) + << " with title = " << std::string( tmphist->GetTitle() ) << std::endl; + if ( tmphist->GetEntries() < 5 ) { + LOGINFO << "Histogram " << std::string(histo) << " has only " << tmphist->GetEntries() + << " entries. Trouble ahead." << std::endl; + } + qbinHist_ [ ii ][ jj ] = tmphist; + qbinGen_ [ ii ][ jj ] = new SimpleHistogramGenerator( tmphist ); + } + + + } + } +} + + + + +//------------------------------------------------------------------------------ +// Destructor. Use file_ pointer to tell whether we loaded the histograms +// from a file (and do not own them), or we built them ourselves and thus need +// to delete them. +//------------------------------------------------------------------------------ +PixelResolutionHistograms::~PixelResolutionHistograms() +{ + //--- Delete histograms, but only if we own them. If + //--- they came from a file, let them be. + // + if ( ! weOwnHistograms_ ) { + //--- Read the histograms from the TFile, the file will take care of them. + file_->Close(); + /// delete file_ ; // no need to delete if unique_ptr<> + /// file_ = 0; + } + else { + //--- We made the histograms, so first write them inthe output ROOT file and close it. + LOGINFO + << "PixelResHistoStore: Writing the histograms to the output file. " // << filename + << std::endl; + file_->Write(); + file_->Close(); + + // ROOT file has the ownership, and once the file is closed, + // all of these guys are deleted. So, we don't need to do anything. + } // else + + + + //--- Delete FastSim generators. (It's safe to delete a nullptr.) + for( int ii=0; ii 2 ? 3 : qbin; + if( icotalpha >= 0 && icotalpha < cotalphaBins_ && icotbeta >= 0 && icotbeta < cotbetaBins_ ) { + qbinHist_[icotbeta][icotalpha]->Fill((double)iqbin); + if( nxpix == 1 ) + resSinglePixelXHist_[icotbeta][icotalpha]->Fill(dx/cmtomicron); + else + resMultiPixelXHist_[icotbeta][icotalpha][iqbin]->Fill(dx/cmtomicron); + if( nypix == 1 ) + resSinglePixelYHist_[icotbeta][icotalpha]->Fill(dy/cmtomicron); + else + resMultiPixelYHist_[icotbeta][icotalpha][iqbin]->Fill(dy/cmtomicron); + } + + return 0; +} + + +//------------------------------------------------------------------------------ +// Return the histogram generator for resolution in X. A generator contains +// both the histogram and knows how to throw a random number off it. It is +// called from FastSim (from PixelTemplateSmearerBase). +// If cotalpha or cotbeta are outside of the range, return the end of the range. +//------------------------------------------------------------------------------ +const SimpleHistogramGenerator * +PixelResolutionHistograms:: +getGeneratorX( double cotalpha, double cotbeta, int qbin, bool single ) +{ + int icotalpha, icotbeta, iqbin ; + icotalpha = (int)floor( (cotalpha - cotalphaLowEdge_) / cotalphaBinWidth_ ) ; + icotbeta = (int)floor( (cotbeta - cotbetaLowEdge_) / cotbetaBinWidth_ ) ; + iqbin = qbin > 2 ? 3 : qbin; // if (qbin>2) then = 3, else return qbin + // + //if( icotalpha >= 0 && icotalpha < cotalphaBins_ && icotbeta >= 0 && icotbeta < cotbetaBins_ ) { + + if (icotalpha < 0) + icotalpha = 0; + if (icotalpha >= cotalphaBins_) + icotalpha = cotalphaBins_ - 1; + + if (icotbeta < 0) + icotbeta = 0; + if (icotbeta >= cotbetaBins_) + icotbeta = cotbetaBins_ - 1; + + // At this point we are sure to return *some bin* from the 3D histogram + + if( single ) + return resSinglePixelXGen_[icotbeta][icotalpha]; + else + return resMultiPixelXGen_[icotbeta][icotalpha][iqbin]; + + // } + //else + //return nullptr; +} + +//------------------------------------------------------------------------------ +// Return the histogram generator for resolution in Y. A generator contains +// both the histogram and knows how to throw a random number off it. It is +// called from FastSim (from PixelTemplateSmearerBase). +// If cotalpha or cotbeta are outside of the range, return the end of the range. +//------------------------------------------------------------------------------ +const SimpleHistogramGenerator * +PixelResolutionHistograms:: +getGeneratorY( double cotalpha, double cotbeta, int qbin, bool single ) +{ + int icotalpha, icotbeta, iqbin ; + icotalpha = (int)floor( (cotalpha - cotalphaLowEdge_) / cotalphaBinWidth_ ) ; + icotbeta = (int)floor( (cotbeta - cotbetaLowEdge_) / cotbetaBinWidth_ ) ; + iqbin = qbin > 2 ? 3 : qbin; // if (qbin>2) then = 3, else return qbin + // + //if( icotalpha >= 0 && icotalpha < cotalphaBins_ && icotbeta >= 0 && icotbeta < cotbetaBins_ ) { + + if (icotalpha < 0) + icotalpha = 0; + if (icotalpha >= cotalphaBins_) + icotalpha = cotalphaBins_ - 1; + + if (icotbeta < 0) + icotbeta = 0; + if (icotbeta >= cotbetaBins_) + icotbeta = cotbetaBins_ - 1; + + // At this point we are sure to return *some bin* from the 3D histogram + + if( single ) + return resSinglePixelYGen_[icotbeta][icotalpha]; + else + return resMultiPixelYGen_[icotbeta][icotalpha][iqbin]; + + //} + //else + //return nullptr; +} diff --git a/FastSimulation/TrackingRecHitProducer/src/PixelTemplateSmearerBase.cc b/FastSimulation/TrackingRecHitProducer/src/PixelTemplateSmearerBase.cc index 7b7e48708d510..7a46c864b840b 100644 --- a/FastSimulation/TrackingRecHitProducer/src/PixelTemplateSmearerBase.cc +++ b/FastSimulation/TrackingRecHitProducer/src/PixelTemplateSmearerBase.cc @@ -13,9 +13,14 @@ #include "FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h" #include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithmFactory.h" #include "FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitProduct.h" +#include "FastSimulation/TrackingRecHitProducer/interface/PixelResolutionHistograms.h" + +// Pixel related stuff +#include "CondFormats/SiPixelObjects/interface/SiPixelTemplateDBObject.h" // Geometry -#include "Geometry/CommonDetUnit/interface/GeomDet.h" +/// #include "Geometry/CommonDetUnit/interface/GeomDetUnit.h" // Keep... needed if we backport to CMSSW_9 +#include "Geometry/TrackerGeometryBuilder/interface/PixelGeomDetUnit.h" #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h" #include "DataFormats/GeometryVector/interface/LocalPoint.h" #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h" @@ -24,13 +29,17 @@ #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h" #include "FastSimulation/Utilities/interface/SimpleHistogramGenerator.h" -// STL +// Framework (includes ESHandle<>) +#include "FWCore/Framework/interface/EventSetup.h" +#include "FWCore/MessageLogger/interface/MessageLogger.h" +#include "CalibTracker/Records/interface/SiPixelTemplateDBObjectESProducerRcd.h" // ROOT #include #include #include +using namespace std; const double microntocm = 0.0001; @@ -40,28 +49,142 @@ PixelTemplateSmearerBase::PixelTemplateSmearerBase( const edm::ParameterSet& config, edm::ConsumesCollector& consumesCollector ): - TrackingRecHitAlgorithm(name,config,consumesCollector) + TrackingRecHitAlgorithm(name,config,consumesCollector) { + //--- Basic stuff mergeHitsOn = config.getParameter("MergeHitsOn"); - templateId = config.getParameter ( "templateId" ); + isBarrel = config.getParameter ( "isBarrel" ); + int detType = (isBarrel) ? 1 : 0; // 1 for barrel, 0 for forward (or could we just promote bool into int...?) + + //--- Resolution file names. + theBigPixelResolutionFileName = config.getParameter( "BigPixelResolutionFile" ); + theBigPixelResolutionFileName = edm::FileInPath( theBigPixelResolutionFileName ).fullPath(); + + theEdgePixelResolutionFileName = config.getParameter( "EdgePixelResolutionFile" ); + theEdgePixelResolutionFileName = edm::FileInPath( theEdgePixelResolutionFileName ).fullPath(); + + theRegularPixelResolutionFileName = config.getParameter( "RegularPixelResolutionFile" ); + theRegularPixelResolutionFileName = edm::FileInPath( theRegularPixelResolutionFileName ).fullPath(); + + + //--- Create the resolution histogram objects, which will load the histograms + // and initialize random number generators. + // + int status = 0; + theRegularPixelResolutions = + std::make_shared( theRegularPixelResolutionFileName, "" ); + if (( status = theRegularPixelResolutions->status()) != 0 ) { + throw cms::Exception("PixelTemplateSmearerBase:") + << " constructing PixelResolutionHistograms file " << theRegularPixelResolutionFileName + << " failed with status = " << status << std::endl; + } + + theBigPixelResolutions = + std::make_shared( theBigPixelResolutionFileName, "", detType, (!isBarrel), false, true ); // can miss qBin + if (( status = theBigPixelResolutions->status()) != 0 ) { + throw cms::Exception("PixelTemplateSmearerBase:") + << " constructing PixelResolutionHistograms file " << theBigPixelResolutionFileName + << " failed with status = " << status << std::endl; + } + + theEdgePixelResolutions = + std::make_shared( theEdgePixelResolutionFileName, "", detType, false, true, true ); // can miss both single & qBin + if (( status = theEdgePixelResolutions->status()) != 0 ) { + throw cms::Exception("PixelTemplateSmearerBase:") + << " constructing PixelResolutionHistograms file " << theEdgePixelResolutionFileName + << " failed with status = " << status << std::endl; + } + + + + //--- Merging info. + theMergingProbabilityFileName = config.getParameter( "MergingProbabilityFile" ); + theMergingProbabilityFileName = edm::FileInPath( theMergingProbabilityFileName ).fullPath(); + theMergingProbabilityFile = std::make_unique( theMergingProbabilityFileName.c_str(), "READ"); + + theMergedPixelResolutionXFileName = config.getParameter( "MergedPixelResolutionXFile" ); + theMergedPixelResolutionXFileName = edm::FileInPath( theMergedPixelResolutionXFileName ).fullPath(); + theMergedPixelResolutionXFile = std::make_unique( theMergedPixelResolutionXFileName.c_str(), "READ"); + + theMergedPixelResolutionYFileName = config.getParameter( "MergedPixelResolutionYFile" ); + theMergedPixelResolutionYFileName = edm::FileInPath( theMergedPixelResolutionYFileName ).fullPath(); + theMergedPixelResolutionYFile = std::make_unique( theMergedPixelResolutionYFileName.c_str(), "READ"); + + + // const SiPixelTemplateDBObject & dbobject; + // const SiPixelTemplateDBObject dbobject; // dummy, just to make it compile &&& + + //--- Load the templates. + if ( config.exists("templateId") ) { + //--- Load template with ID=templateId from a local ascii file. + templateId = config.getParameter ( "templateId" ); + if ( templateId > 0 ) { + if ( !SiPixelTemplate::pushfile(templateId, thePixelTemp_) ) { + throw cms::Exception("PixelTemplateSmearerBase:") + <<"SiPixel Template " << templateId << " Not Loaded Correctly!" << std::endl; + } + } + } + + //--- Else... The templates will be loaded from the DB... + // (They are needed for data and full sim MC, so in a production FastSim + // run, everything should already be in the DB.) + // + // But note that we can do it only at the beginning of the + // event. So nothing happens now. } PixelTemplateSmearerBase::~PixelTemplateSmearerBase() { - for (auto it = theXHistos.begin(); it != theXHistos.end(); ++it ) - { - delete it->second; - } - for (auto it = theYHistos.begin(); it != theYHistos.end(); ++it ) - { - delete it->second; - } - theXHistos.clear(); - theYHistos.clear(); + //--- Delete the templates. This is safe even if thePixelTemp_ vector is empty. + for (auto x : thePixelTemp_) x.destroy(); +} + +//------------------------------------------------------------------------------- +// beginRun(); the templates are loaded in TrackingRecHitProducer, and unpacked +// into the template store. We get their references here, and use them. However, +// if we are loading a dedicated template ID from an ascii file just for this +// rechit smearing algorithm, then we use our own template store. +//------------------------------------------------------------------------------- +void PixelTemplateSmearerBase::beginRun(edm::Run const& run, const edm::EventSetup& eventSetup, + const SiPixelTemplateDBObject * pixelTemplateDBObjectPtr, + std::vector< SiPixelTemplateStore > & tempStoreRef ) +{ + //--- Check if we need to use the template from the DB (namely if + // id == -1). Otherwise the template has already been loaded from + // the ascii file in constructor, and thePixelTempRef wakes up + // pointing to thePixelTemp_, so then we use our own store. + // + if ( templateId == -1 ) { + thePixelTempRef = tempStoreRef; // we use the store from TrackingRecHitProducer + pixelTemplateDBObject_ = pixelTemplateDBObjectPtr; // needed for template<-->DetId map. + } + + + //--- Commented code below (the DB interface) should say here, in case we need it: + // edm::ESHandle templateDBobject; + // eventSetup.get().get(templateDBobject); + // pixelTemplateDBObject_ = templateDBobject.product(); + + // //--- Now that we have the DB object, load the correct templates from the DB. + // // (They are needed for data and full sim MC, so in a production FastSim + // // run, everything should already be in the DB.) + // if ( !SiPixelTemplate::pushfile( *pixelTemplateDBObject_ , thePixelTemp_) ) { + // throw cms::Exception("PixelTemplateSmearerPlugin:") + // <<"SiPixel Template " << templateId << " Not Loaded Correctly!"<smearIt = true; } - // Step 2: iterate over all hits, replace mgbh[j] by mgbh[i] (so that nobody points to i) - MergeGroup * mgbhj = mergeGroupByHit[j]; + // Step 2: iterate over all hits, replace mgbh[j] by mgbh[i] (so that nobody points to i) + MergeGroup * mgbhj = mergeGroupByHit[j]; for ( int k = 0; k < nHits; ++k ) { if ( mgbhj == mergeGroupByHit[k]) { - // Hit k also uses the same merge group, tell them to switch to mgbh[i] + // Hit k also uses the same merge group, tell them to switch to mgbh[i] mergeGroupByHit[k] = mergeGroupByHit[i]; } } mgbhj->smearIt = false; mergeGroupByHit[i]->smearIt = true; - // Step 3 would have been to delete mgbh[j]... however, we'll do that at the end anyway. - // The key was to prevent mgbh[j] from being accessed further, and we have done that, - // since now no mergeGroupByHit[] points to mgbhj any more. Note that the above loop + // Step 3 would have been to delete mgbh[j]... however, we'll do that at the end anyway. + // The key was to prevent mgbh[j] from being accessed further, and we have done that, + // since now no mergeGroupByHit[] points to mgbhj any more. Note that the above loop // also set mergeGroupByHit[i] = mergeGroupByHit[j], too. } } @@ -172,11 +295,9 @@ PixelTemplateSmearerBase::process(TrackingRecHitProductPtr product) const // other hit. Create a new merge group for i and j mergeGroupByHit[i] = new MergeGroup(); listOfMergeGroups.push_back( mergeGroupByHit[i] ); // keep track of it - //std::cout << "new merge group " << mergeGroupByHit[i] << std::endl; // // Add hit i as the first to its own merge group // (simHits[i] is a const pointer to PSimHit). - //std::cout << "ALICE: simHits" << simHits[i] << std::endl; mergeGroupByHit[i]->group.push_back( simHitIdPairs[i] ); mergeGroupByHit[i]->smearIt = true; } @@ -207,9 +328,8 @@ PixelTemplateSmearerBase::process(TrackingRecHitProductPtr product) const } // --- end of if (mergeHitsOn) else { - // std::cout << "Merged Hits Off!" << std::endl; - //Now we've turned off hit merging, so all hits should be pushed - //back to listOfUnmergedHits + // Now we've turned off hit merging, so all hits should be pushed + // back to listOfUnmergedHits for (int i = 0; i < nHits; ++i ) { listOfUnmergedHits.push_back( simHitIdPairs[i] ); @@ -247,7 +367,7 @@ PixelTemplateSmearerBase::process(TrackingRecHitProductPtr product) const //------------------------------------------------------------------------------ // Smear one hit. The main action is in here. //------------------------------------------------------------------------------ -FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit( +FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit ( const PSimHit& simHit, const PixelGeomDetUnit* detUnit, const double boundX, @@ -255,25 +375,32 @@ FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit( RandomEngineAndDistribution const* random) const { - // at the beginning the position is the Local Point in the local pixel module reference frame - // same code as in PixelCPEBase - LocalVector localDir = simHit.momentumAtEntry().unit(); + //--- At the beginning the position is the Local Point in the local pixel module reference frame + // same code as in PixelCPEBase + // + LocalVector localDir = simHit.momentumAtEntry(); // don't need .unit(), we will take the ratio float locx = localDir.x(); float locy = localDir.y(); float locz = localDir.z(); + //--- cotangent of local angles \alpha and \beta. + // alpha: angle with respect to local x axis in local (x,z) plane + // beta: angle with respect to local y axis in local (y,z) plane + // float cotalpha = locx/locz; float cotbeta = locy/locz; - float sign=1.; - - if(isForward) - { - if( cotbeta < 0 ) - { - sign=-1.; - } - cotbeta = sign*cotbeta; - } + + //--- Save the original signs of cot\alpha and cot\beta + int signOfCotalpha = (cotalpha < 0) ? -1 : 1; // sign(cotalpha); + int signOfCotbeta = (cotbeta < 0) ? -1 : 1; // sign(cotbeta); + // + //--- Use absolute values to find the templates from the list + cotalpha *= signOfCotalpha; // = abs(cotalpha) + cotbeta *= signOfCotbeta; // = abs(cotbeta) + + LogDebug ("SmearHit") << "LocalVector=" << locx << "," << locy << "," << locz + << " momentum=" << localDir.mag() + << " cotalpha=" << cotalpha << ", cotbeta=" << cotbeta; const PixelTopology* theSpecificTopology = &(detUnit->specificType().specificTopology()); const RectangularPixelTopology *rectPixelTopology = static_cast(theSpecificTopology); @@ -306,41 +433,54 @@ FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit( int xbin = (int)xhit; float yfrac= yhit - (float)ybin; float xfrac= xhit - (float)xbin; - //Protect againt ybin, xbin being outside of range [0-39] + //Protect againt ybin, xbin being outside of range [0-39] // &&& Why limit of 39? if( ybin < 0 ) ybin = 0; if( ybin > 39 ) ybin = 39; if( xbin < 0 ) xbin = 0; if( xbin > 39 ) xbin = 39; + + int ID = templateId; + if ( templateId == -1 ) { + // We have loaded the whole template set from the DB, + // so ask the DB object to find us the right one. + ID = pixelTemplateDBObject_->getTemplateID( detUnit->geographicalId() ); // need uint32_t detid + // theDetParam.theDet->geographicalId()); + } + + //--- Make the template object + SiPixelTemplate templ( thePixelTempRef ); + + //--- Produce the template that corresponds to our local angles. + templ.interpolate( ID, cotalpha, cotbeta ); + //Variables for SiPixelTemplate output //qBin -- normalized pixel charge deposition float qbin_frac[4]; //Single pixel cluster projection possibility float ny1_frac, ny2_frac, nx1_frac, nx2_frac; bool singlex = false, singley = false; - SiPixelTemplate templ(thePixelTemp_); - templ.interpolate(templateId, cotalpha, cotbeta); - templ.qbin_dist(templateId, cotalpha, cotbeta, qbin_frac, ny1_frac, ny2_frac, nx1_frac, nx2_frac ); + templ.qbin_dist( ID, cotalpha, cotbeta, qbin_frac, ny1_frac, ny2_frac, nx1_frac, nx2_frac ); int nqbin; double xsizeProbability = random->flatShoot(); double ysizeProbability = random->flatShoot(); - bool hitbigx = rectPixelTopology->isItBigPixelInX( (int)mpx ); - bool hitbigy = rectPixelTopology->isItBigPixelInY( (int)mpy ); + bool hitbigx = rectPixelTopology->isItBigPixelInX( (int)mpx ); // pixel we hit in x + bool hitbigy = rectPixelTopology->isItBigPixelInY( (int)mpy ); // pixel we hit in y if( hitbigx ) - if( xsizeProbability < nx2_frac ) singlex = true; - else singlex = false; + if( xsizeProbability < nx2_frac ) singlex = true; + else singlex = false; else - if( xsizeProbability < nx1_frac ) singlex = true; - else singlex = false; + if( xsizeProbability < nx1_frac ) singlex = true; + else singlex = false; if( hitbigy ) - if( ysizeProbability < ny2_frac ) singley = true; - else singley = false; + if( ysizeProbability < ny2_frac ) singley = true; + else singley = false; else - if( ysizeProbability < ny1_frac ) singley = true; - else singley = false; + if( ysizeProbability < ny1_frac ) singley = true; + else singley = false; @@ -383,7 +523,7 @@ FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit( for( firstY = 0; firstY < BYSIZE; ++firstY ) { bool yCluster = ytemp[firstY] > qThreshold; - if(yCluster) + if (yCluster) { offsetY1 = BHY -firstY; break; @@ -392,7 +532,7 @@ FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit( for(lastY = firstY; lastY < BYSIZE; ++lastY) { bool yCluster = ytemp[lastY] > qThreshold; - if(!yCluster) + if (!yCluster) { lastY = lastY - 1; offsetY2 = lastY - BHY; @@ -422,10 +562,10 @@ FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit( //--- Prepare to return results - Local3DPoint thePosition; - double thePositionX; - double thePositionY; - double thePositionZ; + Local3DPoint thePosition; + double theShiftInX; + double theShiftInY; + double theShiftInZ; LocalError theError; double theErrorX; double theErrorY; @@ -459,8 +599,8 @@ FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit( //Variables for SiPixelTemplate pixel hit error output float sigmay, sigmax, sy1, sy2, sx1, sx2; templ.temperrors( - templateId, cotalpha, cotbeta, nqbin, // inputs - sigmay, sigmax, sy1, sy2, sx1, sx2 // outputs + ID, cotalpha, cotbeta, nqbin, // inputs + sigmay, sigmax, sy1, sy2, sx1, sx2 // outputs ); if(edge) @@ -538,115 +678,72 @@ FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit( // Local Error is 2D: (xx,xy,yy), square of sigma in first an third position // as for resolution matrix - // Generate position - // get resolution histograms - int cotalphaHistBin = (int)( ( cotalpha - rescotAlpha_binMin ) / rescotAlpha_binWidth + 1 ); - int cotbetaHistBin = (int)( ( cotbeta - rescotBeta_binMin ) / rescotBeta_binWidth + 1 ); - // protection against out-of-range (undeflows and overflows) - if (cotalphaHistBin < 1) cotalphaHistBin = 1; - if (cotbetaHistBin < 1) cotbetaHistBin = 1; - if (cotalphaHistBin > (int)rescotAlpha_binN) cotalphaHistBin = (int)rescotAlpha_binN; - if (cotbetaHistBin > (int)rescotBeta_binN) cotbetaHistBin = (int)rescotBeta_binN; - // - unsigned int theXHistN; - unsigned int theYHistN; - if (!isForward) - { - if (edge) - { - theXHistN = cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1); - theYHistN = theXHistN; - } - else - { - if (singlex) - { - if (hitbigx) theXHistN = 1 * 100000 + cotalphaHistBin * 100 + cotbetaHistBin ; - else theXHistN = 1 * 10000 + cotbetaHistBin * 10 + cotalphaHistBin ; - } - else - { - if (hasBigPixelInX) theXHistN = 1 * 1000000 + 1 * 100000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1); - else theXHistN = 1 * 100000 + 1 * 10000 + cotbetaHistBin * 100 + cotalphaHistBin * 10 + (nqbin+1); - } - - if(singley) - { - if (hitbigy) theYHistN = 1 * 100000 + cotalphaHistBin * 100 + cotbetaHistBin ; - else theYHistN = 1 * 10000 + cotbetaHistBin * 10 + cotalphaHistBin ; - } - else - { - if (hasBigPixelInY) theYHistN = 1 * 1000000 + 1 * 100000 + cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1); - else theYHistN = 1 * 100000 + 1 * 10000 + cotbetaHistBin * 100 + cotalphaHistBin * 10 + (nqbin+1); - } - } + //--- Next, we need to generate the smeared position. First we need to figure + // out which kind of histograms we are supposed to use for this particular hit. + // These are pointers to the set of histograms used to generate the rec hit + // positions. (We need to handle X and Y separately.) + shared_ptr resHistsX = nullptr; + shared_ptr resHistsY = nullptr; + + + if (edge) { + resHistsX = resHistsY = theEdgePixelResolutions; + singlex = singley = false; // no single resolutions for Edge } - else - { - if (edge) - { - theXHistN = cotalphaHistBin * 1000 + cotbetaHistBin * 10 + (nqbin+1); - theYHistN = theXHistN; - } - else - { - if (singlex) - { - if (hitbigx) theXHistN = 100000 + cotalphaHistBin * 100 + cotbetaHistBin; - else theXHistN = cotbetaHistBin * 10 + cotalphaHistBin; - } - else - { - theXHistN = 10000 + cotbetaHistBin * 100 + cotalphaHistBin * 10 + (nqbin+1); - } - - if(singley) - { - if (hitbigy) theYHistN = 100000 + cotalphaHistBin * 100 + cotbetaHistBin; - else theYHistN = cotbetaHistBin * 10 + cotalphaHistBin; - } - else - { - theYHistN = 10000 + cotbetaHistBin * 100 + cotalphaHistBin * 10 + (nqbin+1); - } - } + else { + //--- Decide resolution histogram set for X + if ( (singlex && hitbigx) || (isBarrel && hasBigPixelInX) ) { + resHistsX = theBigPixelResolutions; + } + else { + resHistsX = theRegularPixelResolutions; + } + //--- Decide resolution histogram set for Y + if ( (singley && hitbigy) || (isBarrel && hasBigPixelInY) ) { + resHistsY = theBigPixelResolutions; + } + else { + resHistsY = theRegularPixelResolutions; + } } + + //--- Get generators, separately for X and for Y. + const SimpleHistogramGenerator * xgen + = resHistsX->getGeneratorX( cotalpha, cotbeta, nqbin, singlex ); + const SimpleHistogramGenerator * ygen + = resHistsY->getGeneratorY( cotalpha, cotbeta, nqbin, singley ); + + //--- Check if we found a histogram. If nullptr, then throw up. + if ( !xgen || !ygen ) { + throw cms::Exception("FastSimulation/TrackingRecHitProducer") + << "Histogram (cot\alpha=" << cotalpha + <<", cot\beta="<< cotbeta << ", nQbin=" << nqbin + << ") was not found for PixelTemplateSmearer. Check if the smearing resolution histogram exists."; + } + - + //--- Smear the hit Position. We do it in the do-while loop in order to + //--- allow multiple tries, in case we generate a rec hit which is outside + //--- of the boundaries of the sensor. unsigned int retry = 0; do { - // - // Smear the hit Position - - std::map::const_iterator xgenIt = theXHistos.find(theXHistN); - std::map::const_iterator ygenIt = theYHistos.find(theYHistN); - if (xgenIt==theXHistos.cend() || ygenIt==theYHistos.cend()) - { - throw cms::Exception("FastSimulation/TrackingRecHitProducer") << "Histogram ("<generate(random); + theShiftInY = ygen->generate(random); - const SimpleHistogramGenerator* xgen = xgenIt->second; - const SimpleHistogramGenerator* ygen = ygenIt->second; + // Now multiply by the sign of the cotangent of appropriate angle + theShiftInX *= signOfCotalpha; + theShiftInY *= signOfCotbeta; - thePositionX = xgen->generate(random); - thePositionY = ygen->generate(random); - - - if( isForward ) - { - thePositionY *= sign; - } - thePositionZ = 0.0; // set at the centre of the active area + theShiftInZ = 0.0; // set to the mid-plane of the sensor. thePosition = Local3DPoint( - simHit.localPosition().x() + thePositionX, - simHit.localPosition().y() + thePositionY, - simHit.localPosition().z() + thePositionZ + simHit.localPosition().x() + theShiftInX, + simHit.localPosition().y() + theShiftInY, + simHit.localPosition().z() + theShiftInZ ); retry++; if (retry > 10) @@ -673,6 +770,11 @@ FastSingleTrackerRecHit PixelTemplateSmearerBase::smearHit( } + + +//------------------------------------------------------------------------------ +// Smear all umerged hits on this DetUnit +//------------------------------------------------------------------------------ TrackingRecHitProductPtr PixelTemplateSmearerBase:: processUnmergedHits( std::vector & unmergedHits, @@ -691,6 +793,10 @@ processUnmergedHits( } + +//------------------------------------------------------------------------------ +// Smear all MERGED hits on this DetUnit +//------------------------------------------------------------------------------ TrackingRecHitProductPtr PixelTemplateSmearerBase:: processMergeGroups( std::vector< MergeGroup* > & mergeGroups, @@ -712,8 +818,12 @@ processMergeGroups( } + +//------------------------------------------------------------------------------ +// Smear all hits MERGED together. This is called a MergeGroup. +//------------------------------------------------------------------------------ FastSingleTrackerRecHit PixelTemplateSmearerBase:: -smearMergeGroup( +smearMergeGroup ( MergeGroup* mg, const PixelGeomDetUnit * detUnit, const double boundX, const double boundY, @@ -750,21 +860,20 @@ smearMergeGroup( float locy = loccy/nHit; float locz = loccz/nHit; - // alpha: angle with respect to local x axis in local (x,z) plane + //--- cotangent of local angles \alpha and \beta. + // alpha: angle with respect to local x axis in local (x,z) plane + // beta: angle with respect to local y axis in local (y,z) plane + // float cotalpha = locx/locz; - // beta: angle with respect to local y axis in local (y,z) plane float cotbeta = locy/locz; - float sign=1.; - - if( isForward ) - { - if( cotbeta < 0 ) - { - sign=-1.; - } - cotbeta = sign*cotbeta; - } + //--- Save the original signs of cot\alpha and cot\beta + int signOfCotalpha = (cotalpha < 0) ? -1 : 1; // sign(cotalpha); + int signOfCotbeta = (cotbeta < 0) ? -1 : 1; // sign(cotbeta); + // + //--- Use absolute values to find the templates from the list + cotalpha *= signOfCotalpha; // = abs(cotalpha) + cotbeta *= signOfCotbeta; // = abs(cotbeta) float lpx = locpx/nHit; float lpy = locpy/nHit; @@ -782,21 +891,33 @@ smearMergeGroup( int xbin = (int)xhit; float yfrac= yhit - (float)ybin; float xfrac= xhit - (float)xbin; - //Protect againt ybin, xbin being outside of range [0-39] + // Protect againt ybin, xbin being outside of range [0-39] if( ybin < 0 ) ybin = 0; if( ybin > 39 ) ybin = 39; if( xbin < 0 ) xbin = 0; if( xbin > 39 ) xbin = 39; - //Variables for SiPixelTemplate output - //qBin -- normalized pixel charge deposition + int ID = templateId; + if ( templateId == -1 ) { + // We have loaded the whole template set from the DB, + // so ask the DB object to find us the right one. + ID = pixelTemplateDBObject_->getTemplateID( detUnit->geographicalId() ); // need uint32_t detid + // theDetParam.theDet->geographicalId()); + } + + //--- Make the template object + SiPixelTemplate templ( thePixelTempRef ); + + //--- Produce the template that corresponds to our local angles. + templ.interpolate( ID, cotalpha, cotbeta ); + + // Variables for SiPixelTemplate output + // qBin -- normalized pixel charge deposition float qbin_frac[4]; - //Single pixel cluster projection possibility + // Single pixel cluster projection possibility float ny1_frac, ny2_frac, nx1_frac, nx2_frac; bool singlex = false, singley = false; - SiPixelTemplate templ(thePixelTemp_); - templ.interpolate(templateId, cotalpha, cotbeta); - templ.qbin_dist(templateId, cotalpha, cotbeta, qbin_frac, ny1_frac, ny2_frac, nx1_frac, nx2_frac ); + templ.qbin_dist( ID, cotalpha, cotbeta, qbin_frac, ny1_frac, ny2_frac, nx1_frac, nx2_frac ); int nqbin; // double xsizeProbability = random->flatShoot(); @@ -807,6 +928,7 @@ smearMergeGroup( // random multiplicity for alpha and beta + double qbinProbability = random->flatShoot(); for(int i = 0; i<4; ++i) { @@ -834,14 +956,12 @@ smearMergeGroup( //--- Prepare to return results Local3DPoint thePosition; - double thePositionX; - double thePositionY; - double thePositionZ; + double theShiftInX; + double theShiftInY; + double theShiftInZ; LocalError theError; double theErrorX; double theErrorY; - //double theErrorZ; - //------------------------------ @@ -856,7 +976,7 @@ smearMergeGroup( //Variables for SiPixelTemplate pixel hit error output float sigmay, sigmax, sy1, sy2, sx1, sx2; - templ.temperrors(templateId, cotalpha, cotbeta, nqbin, // inputs + templ.temperrors( ID, cotalpha, cotbeta, nqbin, // inputs sigmay, sigmax, sy1, sy2, sx1, sx2 ); // outputs // define private mebers --> Errors @@ -923,18 +1043,20 @@ smearMergeGroup( const SimpleHistogramGenerator* xgen = new SimpleHistogramGenerator( (TH1F*) theMergedPixelResolutionXFile-> Get("th1x")); const SimpleHistogramGenerator* ygen = new SimpleHistogramGenerator( (TH1F*) theMergedPixelResolutionYFile-> Get("th1y")); - thePositionX = xgen->generate(random); - thePositionY = ygen->generate(random); + // Generate the position (x,y of the rec hit). + theShiftInX = xgen->generate(random); + theShiftInY = ygen->generate(random); - if( isForward ) - { - thePositionY *= sign; - } - thePositionZ = 0.0; // set at the centre of the active area - thePosition = - Local3DPoint(lpx + thePositionX , - lpy + thePositionY , - lpz + thePositionZ ); + // Now multiply by the sign of the cotangent of appropriate angle + theShiftInX *= signOfCotalpha; + theShiftInY *= signOfCotbeta; + + theShiftInZ = 0.0; // set at the centre of the active area + + thePosition = + Local3DPoint( lpx + theShiftInX , + lpy + theShiftInY , + lpz + theShiftInZ ); retry++; if (retry > 10) @@ -977,23 +1099,3 @@ bool PixelTemplateSmearerBase::hitsMerge(const PSimHit& simHit1,const PSimHit& s } -//----------------------------------------------------------------------------- -// The isFlipped() is a silly way to determine which detectors are inverted. -// In the barrel for every 2nd ladder the E field direction is in the -// global r direction (points outside from the z axis), every other -// ladder has the E field inside. Something similar is in the -// forward disks (2 sides of the blade). This has to be recognised -// because the charge sharing effect is different. -// -// The isFliped does it by looking and the relation of the local (z always -// in the E direction) to global coordinates. There is probably a much -// better way.(PJ: And faster!) -//----------------------------------------------------------------------------- -bool PixelTemplateSmearerBase::isFlipped(const PixelGeomDetUnit* theDet) const -{ - float tmp1 = theDet->surface().toGlobal(Local3DPoint(0.,0.,0.)).perp(); - float tmp2 = theDet->surface().toGlobal(Local3DPoint(0.,0.,1.)).perp(); - return tmp2(id); } + +void TrackingRecHitAlgorithm::beginRun(edm::Run const& run, const edm::EventSetup& eventSetup, + const SiPixelTemplateDBObject * pixelTemplateDBObjectPtr, + std::vector< SiPixelTemplateStore > & tempStoreRef ) +{ + // The default is to do nothing. +} + + void TrackingRecHitAlgorithm::beginEvent(edm::Event& event, const edm::EventSetup& eventSetup) { edm::ESHandle trackerTopologyHandle; From 34931aeb0f02c05ac2929de4bacbf9311c19281b Mon Sep 17 00:00:00 2001 From: Petar Date: Sat, 29 Sep 2018 20:28:52 -0500 Subject: [PATCH 2/2] Synchronizing rebase of 10.3.0 and 9.4.X --- .../TrackingRecHitProducer/interface/TrackingRecHitAlgorithm.h | 2 +- .../TrackingRecHitProducer/plugins/TrackingRecHitProducer.cc | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithm.h b/FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithm.h index 9c381028c73fe..6d1f24ca20fef 100644 --- a/FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithm.h +++ b/FastSimulation/TrackingRecHitProducer/interface/TrackingRecHitAlgorithm.h @@ -10,7 +10,7 @@ // Pixel-related stuff: #include "CondFormats/SiPixelObjects/interface/SiPixelTemplateDBObject.h" -#include "CondFormats/SiPixelTransient/interface/SiPixelTemplate.h" +#include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplate.h" #include #include diff --git a/FastSimulation/TrackingRecHitProducer/plugins/TrackingRecHitProducer.cc b/FastSimulation/TrackingRecHitProducer/plugins/TrackingRecHitProducer.cc index be9c780161648..f512196eb3962 100644 --- a/FastSimulation/TrackingRecHitProducer/plugins/TrackingRecHitProducer.cc +++ b/FastSimulation/TrackingRecHitProducer/plugins/TrackingRecHitProducer.cc @@ -30,7 +30,7 @@ #include "Geometry/Records/interface/TrackerTopologyRcd.h" // Pixel-related stuff: #include "CondFormats/SiPixelObjects/interface/SiPixelTemplateDBObject.h" -#include "CondFormats/SiPixelTransient/interface/SiPixelTemplate.h" +#include "RecoLocalTracker/SiPixelRecHits/interface/SiPixelTemplate.h" #include "CalibTracker/Records/interface/SiPixelTemplateDBObjectESProducerRcd.h" #include "FastSimulation/TrackingRecHitProducer/interface/TrackerDetIdSelector.h"