New pixel RecHit infrastructure for FastSim (for Phases 0, 1, 2, and beyond)

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 <memory>
+
+/// #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( const char * filename,    // ROOT file for histograms
+			     const char * rootdir,     // Subdirectory in the file, "" if none
+			     const char * 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( const char * filename,      // ROOT file for histograms
+			     const char * 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<TFile> 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

FastSimulation/TrackingRecHitProducer/interface/PixelTemplateSmearerBase.h

@@ -27,14 +27,15 @@
 #include "DataFormats/GeometryVector/interface/Point3DBase.h"
 #include "DataFormats/GeometrySurface/interface/LocalError.h"
 
-// STL
+// STL.  <memory> needed for uniq_ptr<>
 #include <vector>
 #include <string>
 #include <memory>
 
 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<unsigned int, const SimpleHistogramGenerator*> theXHistos;
-        std::map<unsigned int, const SimpleHistogramGenerator*> theYHistos;
-
-        std::unique_ptr<TFile> theEdgePixelResolutionFile;
+	//--- The histogram storage containers.
+        std::shared_ptr<PixelResolutionHistograms> theEdgePixelResolutions;
         std::string theEdgePixelResolutionFileName;
-        std::unique_ptr<TFile> theBigPixelResolutionFile;
+
+        std::shared_ptr<PixelResolutionHistograms> theBigPixelResolutions;
         std::string theBigPixelResolutionFileName;
-        std::unique_ptr<TFile> theRegularPixelResolutionFile;
+
+        std::shared_ptr<PixelResolutionHistograms> theRegularPixelResolutions;
         std::string theRegularPixelResolutionFileName;
+
+	//--- Files with hit merging information:
         std::unique_ptr<TFile> theMergingProbabilityFile;
         std::string theMergingProbabilityFileName;
+
         std::unique_ptr<TFile> theMergedPixelResolutionXFile;
         std::string theMergedPixelResolutionXFileName;
-        std::unique_ptr<TFile> theMergedPixelResolutionYFile;                                                                                        
-        std::string theMergedPixelResolutionYFileName;
 
-        unsigned int theLayer;
+        std::unique_ptr<TFile> 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( 

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 <string>
 #include <memory>
 
@@ -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