TrackingRecHitProducer/src/PixelResolutionHistograms.cc

0001 //
0002 //  The implementation of the PixelResolutinHistograms.cc class.  Please
0003 //  look at PixelResolutionHistograms.h header file for the interface.
0004 //
0005 //------------------------------------------------------------------------------
0006
0007 // The switch, undefined in CMSSW release, and defined by standalone compilation script:
0008
0009 #ifdef SI_PIXEL_TEMPLATE_STANDALONE
0010 //
0011 //--- Stand-alone: Include a the header file from the local directory, as well as
0012 //    dummy implementations of SimpleHistogramGenerator, LogInfo, LogError and LogDebug...
0013 //
0014 #include "PixelResolutionHistograms.h"
0015 //
0016 class TH1F;
0017 class TH2F;
0018 class SimpleHistogramGenerator {
0019 public:
0020   SimpleHistogramGenerator(TH1F* hist) : hist_(hist){};
0021
0022 private:
0023   TH1F* hist_;  // we don't own it
0024 };
0025 #define LOGDEBUG std::cout
0026 #define LOGERROR std::cout
0027 #define LOGINFO std::cout
0028 //
0029 #else
0030 //--- We're inside a CMSSW release: Include the real thing.
0031 //
0032 #include "FastSimulation/TrackingRecHitProducer/interface/PixelResolutionHistograms.h"
0033 #include "FastSimulation/Utilities/interface/RandomEngineAndDistribution.h"
0034 #include "FastSimulation/Utilities/interface/SimpleHistogramGenerator.h"
0035 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0036 //
0037 #define LOGDEBUG LogDebug("")
0038 #define LOGERROR edm::LogError("Error")
0039 #define LOGINFO edm::LogInfo("Info")
0040 //
0041 #endif
0042
0043 // Generic C stuff
0044 #include <cmath>
0045 #include <iostream>
0046 #include <string>
0047
0048 // ROOT
0049 #include <TFile.h>
0050 #include <TH1F.h>
0051 #include <TH2F.h>
0052
0053 // Global definitions
0054 const float cmtomicron = 10000.0;
0055
0056 //------------------------------------------------------------------------------
0057 //  Constructor: Books the FastSim Histograms, given the input parameters
0058 //  which are provided as arguments. These variables are then const inside
0059 //  the class. (That is, once we make the histograms, we can't change the
0060 //  definition of the binning.)
0061 //------------------------------------------------------------------------------
0062 PixelResolutionHistograms::PixelResolutionHistograms(std::string filename,   // ROOT file for histograms
0063                                                      std::string rootdir,    // Subdirectory in the file, "" if none
0064                                                      std::string descTitle,  // Descriptive title
0065                                                      unsigned int detType,   // Where we are... (&&& do we need this?)
0066                                                      double cotbetaBinWidth,
0067                                                      double cotbetaLowEdge,
0068                                                      int cotbetaBins,
0069                                                      double cotalphaBinWidth,
0070                                                      double cotalphaLowEdge,
0071                                                      int cotalphaBins)  //,
0072                                                                         //int   qbinWidth,
0073                                                                         //int   qbins )
0074     : weOwnHistograms_(true),                                           // we'll be making some histos
0075       detType_(detType),
0076       cotbetaBinWidth_(cotbetaBinWidth),
0077       cotbetaLowEdge_(cotbetaLowEdge),
0078       cotbetaBins_(cotbetaBins),
0079       cotalphaBinWidth_(cotalphaBinWidth),
0080       cotalphaLowEdge_(cotalphaLowEdge),
0081       cotalphaBins_(cotalphaBins),
0082       qbinWidth_(1),
0083       qbins_(4),
0084       binningHisto_(nullptr),
0085       resMultiPixelXHist_(),
0086       resSinglePixelXHist_(),  // all to nullptr
0087       resMultiPixelYHist_(),
0088       resSinglePixelYHist_(),  // all to nullptr
0089       qbinHist_(),             // all to nullptr
0090       file_(nullptr),
0091       status_(0),
0092       resMultiPixelXGen_(),
0093       resSinglePixelXGen_(),
0094       resMultiPixelYGen_(),
0095       resSinglePixelYGen_(),
0096       qbinGen_() {
0097   file_ = std::make_unique<TFile>(filename.c_str(), "RECREATE");
0098   //Resolution binning
0099   // const double   cotbetaBinWidth = 1.0;
0100   // const double   cotbetaLowEdge  = -11.5 ;
0101   // const int  cotbetaBins = 23;
0102   // const double   cotalphaBinWidth    = 0.08 ;
0103   // const double   cotalphaLowEdge = -0.36 ;
0104   // const int  cotalphaBins    = 9;
0105   // const int  qbinWidth   = 1;
0106   // const int  qbins       = 4;
0107
0108   // Dummy 2D histogram to store binning:
0109   binningHisto_ = new TH2F("ResHistoBinning",
0110                            descTitle.c_str(),
0111                            cotbetaBins,
0112                            cotbetaLowEdge,
0113                            cotbetaLowEdge + cotbetaBins * cotbetaBinWidth,
0114                            cotalphaBins,
0115                            cotalphaLowEdge,
0116                            cotalphaLowEdge + cotalphaBins * cotalphaBinWidth);
0117
0118   // Store detType in the underflow bin
0119   binningHisto_->SetBinContent(0, 0, detType_);
0120
0121   // All other histograms:
0122   Char_t histo[200];
0123   Char_t title[200];
0124   //
0125   //--- Histograms for clusters with multiple pixels hit in a given direction.
0126   //
0127   for (int ii = 0; ii < cotbetaBins_; ii++) {
0128     for (int jj = 0; jj < cotalphaBins_; jj++) {
0129       for (int kk = 0; kk < qbins_; kk++) {
0130         //
0131         sprintf(histo, "hx%d1%02d%d%d", detType_, ii + 1, jj + 1, kk + 1);  //information of bits of histogram names
0132         //--- First bit 1/0 barrel/forward, second 1/0 multi/single, cotbeta, cotalpha, qbins
0133         sprintf(title,
0134                 "cotbeta %.1f-%.1f cotalpha %.2f-%.2f qbin %d npixel>1 X",
0135                 cotbetaLowEdge_ + ii * cotbetaBinWidth_,
0136                 cotbetaLowEdge_ + (ii + 1) * cotbetaBinWidth_,
0137                 cotalphaLowEdge_ + jj * cotalphaBinWidth_,
0138                 cotalphaLowEdge_ + (jj + 1) * cotalphaBinWidth_,
0139                 kk + 1);
0140         //
0141         resMultiPixelXHist_[ii][jj][kk] = new TH1F(histo, title, 1000, -0.05, 0.05);
0142
0143         sprintf(histo, "hy%d1%02d%d%d", detType_, ii + 1, jj + 1, kk + 1);
0144         sprintf(title,
0145                 "cotbeta %.1f-%.1f cotalpha %.2f-%.2f qbin %d npixel>1 Y",
0146                 cotbetaLowEdge_ + ii * cotbetaBinWidth_,
0147                 cotbetaLowEdge_ + (ii + 1) * cotbetaBinWidth_,
0148                 cotalphaLowEdge_ + jj * cotalphaBinWidth_,
0149                 cotalphaLowEdge_ + (jj + 1) * cotalphaBinWidth_,
0150                 kk + 1);
0151         //
0152         resMultiPixelYHist_[ii][jj][kk] = new TH1F(histo, title, 1000, -0.05, 0.05);
0153       }
0154     }
0155   }
0156
0157   //
0158   //--- Histograms for clusters where only a single pixel was hit in a given direction.
0159   //
0160   for (int ii = 0; ii < cotbetaBins_; ii++) {
0161     for (int jj = 0; jj < cotalphaBins_; jj++) {
0162       sprintf(histo, "hx%d0%02d%d", detType_, ii + 1, jj + 1);  //information of bits of histogram names
0163       //first bit 1/0 barrel/forward, second 1/0 multi/single, cotbeta, cotalpha
0164       sprintf(title,
0165               "cotbeta %.1f-%.1f cotalpha %.2f-%.2f npixel=1 X",
0166               cotbetaLowEdge_ + ii * cotbetaBinWidth_,
0167               cotbetaLowEdge_ + (ii + 1) * cotbetaBinWidth_,
0168               cotalphaLowEdge_ + jj * cotalphaBinWidth_,
0169               cotalphaLowEdge_ + (jj + 1) * cotalphaBinWidth_);
0170       //
0171       resSinglePixelXHist_[ii][jj] = new TH1F(histo, title, 1000, -0.05, 0.05);
0172
0173       sprintf(histo, "hy%d0%02d%d", detType_, ii + 1, jj + 1);
0174       sprintf(title,
0175               "cotbeta %.1f-%.1f cotalpha %.2f-%.2f npixel=1 Y",
0176               cotbetaLowEdge_ + ii * cotbetaBinWidth_,
0177               cotbetaLowEdge_ + (ii + 1) * cotbetaBinWidth_,
0178               cotalphaLowEdge_ + jj * cotalphaBinWidth_,
0179               cotalphaLowEdge_ + (jj + 1) * cotalphaBinWidth_);
0180       //
0181       resSinglePixelYHist_[ii][jj] = new TH1F(histo, title, 1000, -0.05, 0.05);
0182
0183       sprintf(histo, "hqbin%d%02d%d", detType_, ii + 1, jj + 1);
0184       sprintf(title,
0185               "cotbeta %.1f-%.1f cotalpha %.2f-%.2f qbin",
0186               cotbetaLowEdge_ + ii * cotbetaBinWidth_,
0187               cotbetaLowEdge_ + (ii + 1) * cotbetaBinWidth_,
0188               cotalphaLowEdge_ + jj * cotalphaBinWidth_,
0189               cotalphaLowEdge_ + (jj + 1) * cotalphaBinWidth_);
0190       //
0191       qbinHist_[ii][jj] = new TH1F(histo, title, 4, -0.49, 3.51);
0192     }
0193   }
0194 }
0195
0196 //------------------------------------------------------------------------------
0197 //  Another constructor: load the histograms from one file.
0198 //     filename = full path to filename
0199 //     rootdir  = ROOT directory inside the file
0200 //
0201 //  The other parameters are the same (needed later) and must correspond
0202 //  to the histograms we are loading from the file.
0203 //------------------------------------------------------------------------------
0204 PixelResolutionHistograms::PixelResolutionHistograms(std::string filename,
0205                                                      std::string rootdir,
0206                                                      int detType,
0207                                                      bool ignore_multi,
0208                                                      bool ignore_single,
0209                                                      bool ignore_qBin)
0210     : weOwnHistograms_(false),  // resolution histograms are owned by the ROOT file
0211       detType_(-1),
0212       cotbetaBinWidth_(0),
0213       cotbetaLowEdge_(0),
0214       cotbetaBins_(0),
0215       cotalphaBinWidth_(0),
0216       cotalphaLowEdge_(0),
0217       cotalphaBins_(0),
0218       qbinWidth_(1),
0219       qbins_(4),
0220       binningHisto_(nullptr),
0221       resMultiPixelXHist_(),
0222       resSinglePixelXHist_(),  // all to nullptr
0223       resMultiPixelYHist_(),
0224       resSinglePixelYHist_(),  // all to nullptr
0225       qbinHist_(),             // all to nullptr
0226       file_(nullptr),
0227       status_(0),
0228       resMultiPixelXGen_(),
0229       resSinglePixelXGen_(),
0230       resMultiPixelYGen_(),
0231       resSinglePixelYGen_(),
0232       qbinGen_() {
0233   Char_t histo[200];        // the name of the histogram
0234   Char_t title[200];        // histo title, for debugging and sanity checking (compare inside file)
0235   TH1F* tmphist = nullptr;  // cache for histo pointer
0236
0237   //--- Open the file for reading.
0238   file_ = std::make_unique<TFile>(filename.c_str(), "READ");
0239   if (!file_) {
0240     status_ = 1;
0241     LOGERROR << "PixelResolutionHistograms:: Error, file " << filename << " not found.";
0242     return;  // PixelTemplateSmearerBase will throw an exception upon our return.
0243   }
0244
0245   //--- The dummy 2D histogram with the binning of cot\beta and cot\alpha:
0246   binningHisto_ = (TH2F*)file_->Get(Form("%s%s", rootdir.c_str(), "ResHistoBinning"));
0247   if (!binningHisto_) {
0248     status_ = 11;
0249     LOGERROR << "PixelResolutionHistograms:: Error, binning histogrram ResHistoBinning not found.";
0250     return;  // PixelTemplateSmearerBase will throw an exception upon our return.
0251   }
0252
0253   if (detType == -1) {
0254     //--- Fish out detType from the underflow bin:
0255     detType_ = binningHisto_->GetBinContent(0, 0);
0256   } else {
0257     detType_ = detType;  // constructor's argument overrides what's in ResHistoBinning histogram.
0258   }
0259
0260   //--- Now we fill the binning variables:
0261   cotbetaAxis_ = binningHisto_->GetXaxis();
0262   cotbetaBinWidth_ = binningHisto_->GetXaxis()->GetBinWidth(1);  // assume all same width
0263   cotbetaLowEdge_ = binningHisto_->GetXaxis()->GetXmin();        // low edge of the first bin
0264   cotbetaBins_ = binningHisto_->GetXaxis()->GetNbins();
0265   cotalphaAxis_ = binningHisto_->GetYaxis();
0266   cotalphaBinWidth_ = binningHisto_->GetYaxis()->GetBinWidth(1);  // assume all same width;
0267   cotalphaLowEdge_ = binningHisto_->GetYaxis()->GetXmin();        // low edge of the first bin;
0268   cotalphaBins_ = binningHisto_->GetYaxis()->GetNbins();
0269
0270   if (!ignore_multi) {
0271     //
0272     //--- Histograms for clusters with multiple pixels hit in a given direction.
0273     //
0274     for (int ii = 0; ii < cotbetaBins_; ii++) {
0275       for (int jj = 0; jj < cotalphaBins_; jj++) {
0276         for (int kk = 0; kk < qbins_; kk++) {
0277           //
0278           sprintf(histo, "hx%d1%02d%d%d", detType_, ii + 1, jj + 1, kk + 1);  //information of bits of histogram names
0279           //--- First bit 1/0 barrel/forward, second 1/0 multi/single, cotbeta, cotalpha, qbins
0280           sprintf(title,
0281                   "cotbeta %.1f-%.1f cotalpha %.2f-%.2f qbin %d npixel>1 X",
0282                   cotbetaLowEdge_ + ii * cotbetaBinWidth_,
0283                   cotbetaLowEdge_ + (ii + 1) * cotbetaBinWidth_,
0284                   cotalphaLowEdge_ + jj * cotalphaBinWidth_,
0285                   cotalphaLowEdge_ + (jj + 1) * cotalphaBinWidth_,
0286                   kk + 1);
0287           //
0288           tmphist = (TH1F*)file_->Get(Form("%s%s", rootdir.c_str(), histo));
0289           if (!tmphist) {
0290             status_ = 2;
0291             LOGERROR << "Failed to find histogram=" << std::string(histo);
0292             return;
0293           }
0294           LOGDEBUG << "Found histo " << std::string(histo) << " with title = " << std::string(tmphist->GetTitle())
0295                    << std::endl;
0296           if (tmphist->GetEntries() < 5) {
0297             LOGINFO << "Histogram " << std::string(histo) << " has only " << tmphist->GetEntries()
0298                     << " entries. Trouble ahead." << std::endl;
0299           }
0300           resMultiPixelXHist_[ii][jj][kk] = tmphist;
0301           resMultiPixelXGen_[ii][jj][kk] = new SimpleHistogramGenerator(tmphist);
0302
0303           sprintf(histo, "hy%d1%02d%d%d", detType_, ii + 1, jj + 1, kk + 1);
0304           sprintf(title,
0305                   "cotbeta %.1f-%.1f cotalpha %.2f-%.2f qbin %d npixel>1 Y",
0306                   cotbetaLowEdge_ + ii * cotbetaBinWidth_,
0307                   cotbetaLowEdge_ + (ii + 1) * cotbetaBinWidth_,
0308                   cotalphaLowEdge_ + jj * cotalphaBinWidth_,
0309                   cotalphaLowEdge_ + (jj + 1) * cotalphaBinWidth_,
0310                   kk + 1);
0311           //
0312           tmphist = (TH1F*)file_->Get(Form("%s%s", rootdir.c_str(), histo));
0313           if (!tmphist) {
0314             status_ = 3;
0315             LOGERROR << "Failed to find histogram=" << std::string(histo);
0316             return;
0317           }
0318           LOGDEBUG << "Found histo " << std::string(histo) << " with title = " << std::string(tmphist->GetTitle())
0319                    << std::endl;
0320           if (tmphist->GetEntries() < 5) {
0321             LOGINFO << "Histogram " << std::string(histo) << " has only " << tmphist->GetEntries()
0322                     << " entries. Trouble ahead." << std::endl;
0323           }
0324           resMultiPixelYHist_[ii][jj][kk] = tmphist;
0325           resMultiPixelYGen_[ii][jj][kk] = new SimpleHistogramGenerator(tmphist);
0326         }
0327       }
0328     }
0329     //
0330   }  // if (not ignore multi)
0331
0332   //
0333   //--- Histograms for clusters where only a single pixel was hit in a given direction.
0334   //
0335   for (int ii = 0; ii < cotbetaBins_; ii++) {
0336     for (int jj = 0; jj < cotalphaBins_; jj++) {
0337       //--- Single pixel, along X.
0338       //
0339       sprintf(histo, "hx%d0%02d%d", detType_, ii + 1, jj + 1);  //information of bits of histogram names
0340       //--- First bit 1/0 barrel/forward, second 1/0 multi/single, cotbeta, cotalpha
0341       sprintf(title,
0342               "cotbeta %.1f-%.1f cotalpha %.2f-%.2f npixel=1 X",
0343               cotbetaLowEdge_ + ii * cotbetaBinWidth_,
0344               cotbetaLowEdge_ + (ii + 1) * cotbetaBinWidth_,
0345               cotalphaLowEdge_ + jj * cotalphaBinWidth_,
0346               cotalphaLowEdge_ + (jj + 1) * cotalphaBinWidth_);
0347       //
0348       tmphist = (TH1F*)file_->Get(Form("%s%s", rootdir.c_str(), histo));
0349       if (!tmphist) {
0350         if (!ignore_single) {
0351           LOGERROR << "Failed to find histogram=" << std::string(histo);
0352           status_ = 4;
0353           return;
0354         }
0355       } else {
0356         LOGDEBUG << "Found histo " << std::string(histo) << " with title = " << std::string(tmphist->GetTitle())
0357                  << std::endl;
0358         LOGDEBUG << "Found histo with title = " << std::string(tmphist->GetTitle()) << std::endl;
0359         if (tmphist->GetEntries() < 5) {
0360           LOGINFO << "Histogram " << std::string(histo) << " has only " << tmphist->GetEntries()
0361                   << " entries. Trouble ahead." << std::endl;
0362         }
0363         resSinglePixelXHist_[ii][jj] = tmphist;
0364         resSinglePixelXGen_[ii][jj] = new SimpleHistogramGenerator(tmphist);
0365       }
0366
0367       //--- Single pixel, along Y.
0368       //
0369       sprintf(histo, "hy%d0%02d%d", detType_, ii + 1, jj + 1);
0370       sprintf(title,
0371               "cotbeta %.1f-%.1f cotalpha %.2f-%.2f npixel=1 Y",
0372               cotbetaLowEdge_ + ii * cotbetaBinWidth_,
0373               cotbetaLowEdge_ + (ii + 1) * cotbetaBinWidth_,
0374               cotalphaLowEdge_ + jj * cotalphaBinWidth_,
0375               cotalphaLowEdge_ + (jj + 1) * cotalphaBinWidth_);
0376       //
0377       tmphist = (TH1F*)file_->Get(Form("%s%s", rootdir.c_str(), histo));
0378       if (!tmphist) {
0379         if (!ignore_single) {
0380           LOGERROR << "Failed to find histogram=" << std::string(histo);
0381           status_ = 5;
0382           return;
0383         }
0384       } else {
0385         LOGDEBUG << "Found histo " << std::string(histo) << " with title = " << std::string(tmphist->GetTitle())
0386                  << std::endl;
0387         if (tmphist->GetEntries() < 5) {
0388           LOGINFO << "Histogram " << std::string(histo) << " has only " << tmphist->GetEntries()
0389                   << " entries. Trouble ahead." << std::endl;
0390         }
0391         resSinglePixelYHist_[ii][jj] = tmphist;
0392         resSinglePixelYGen_[ii][jj] = new SimpleHistogramGenerator(tmphist);
0393       }
0394
0395       //--- qBin distribution, for this (cotbeta, cotalpha) bin.
0396       //
0397       sprintf(histo, "hqbin%d%02d%d", detType_, ii + 1, jj + 1);
0398       sprintf(title,
0399               "cotbeta %.1f-%.1f cotalpha %.2f-%.2f qbin",
0400               cotbetaLowEdge_ + ii * cotbetaBinWidth_,
0401               cotbetaLowEdge_ + (ii + 1) * cotbetaBinWidth_,
0402               cotalphaLowEdge_ + jj * cotalphaBinWidth_,
0403               cotalphaLowEdge_ + (jj + 1) * cotalphaBinWidth_);
0404       //
0405       tmphist = (TH1F*)file_->Get(Form("%s%s", rootdir.c_str(), histo));
0406       if (!tmphist) {
0407         if (!ignore_qBin) {
0408           LOGERROR << "Failed to find histogram=" << std::string(histo);
0409           status_ = 6;
0410           return;
0411         }
0412       } else {
0413         LOGDEBUG << "Found histo " << std::string(histo) << " with title = " << std::string(tmphist->GetTitle())
0414                  << std::endl;
0415         if (tmphist->GetEntries() < 5) {
0416           LOGINFO << "Histogram " << std::string(histo) << " has only " << tmphist->GetEntries()
0417                   << " entries. Trouble ahead." << std::endl;
0418         }
0419         qbinHist_[ii][jj] = tmphist;
0420         qbinGen_[ii][jj] = new SimpleHistogramGenerator(tmphist);
0421       }
0422     }
0423   }
0424 }
0425
0426 //------------------------------------------------------------------------------
0427 //  Destructor.  Use file_ pointer to tell whether we loaded the histograms
0428 //  from a file (and do not own them), or we built them ourselves and thus need
0429 //  to delete them.
0430 //------------------------------------------------------------------------------
0431 PixelResolutionHistograms::~PixelResolutionHistograms() {
0432   //--- Delete histograms, but only if we own them. If
0433   //--- they came from a file, let them be.
0434   //
0435   if (!weOwnHistograms_) {
0436     //--- Read the histograms from the TFile, the file will take care of them.
0437     file_->Close();
0438     /// delete file_ ;   // no need to delete if unique_ptr<>
0439     /// file_ = 0;
0440   } else {
0441     //--- We made the histograms, so first write them inthe output ROOT file and close it.
0442     LOGINFO << "PixelResHistoStore: Writing the histograms to the output file. "  // << filename
0443             << std::endl;
0444     file_->Write();
0445     file_->Close();
0446
0447     // ROOT file has the ownership, and once the file is closed,
0448     // all of these guys are deleted.  So, we don't need to do anything.
0449   }  // else
0450
0451   //--- Delete FastSim generators. (It's safe to delete a nullptr.)
0452   for (int ii = 0; ii < cotbetaBins_; ii++) {
0453     for (int jj = 0; jj < cotalphaBins_; jj++) {
0454       for (int kk = 0; kk < qbins_; kk++) {
0455         delete resMultiPixelXGen_[ii][jj][kk];
0456         delete resMultiPixelYGen_[ii][jj][kk];
0457       }
0458     }
0459   }
0460   for (int ii = 0; ii < cotbetaBins_; ii++) {
0461     for (int jj = 0; jj < cotalphaBins_; jj++) {
0462       delete resSinglePixelXGen_[ii][jj];
0463       delete resSinglePixelYGen_[ii][jj];
0464       delete qbinGen_[ii][jj];
0465     }
0466   }
0467 }
0468
0469 //------------------------------------------------------------------------------
0470 //  Fills the appropriate FastSim histograms.
0471 //  Returns 0 if the relevant histogram(s) were found and filled, 1 if not.
0472 //------------------------------------------------------------------------------
0473 int PixelResolutionHistograms::Fill(
0474     double dx, double dy, double cotalpha, double cotbeta, int qbin, int nxpix, int nypix) {
0475   int icotalpha, icotbeta, iqbin;
0476   icotalpha = (int)floor((cotalpha - cotalphaLowEdge_) / cotalphaBinWidth_);
0477   icotbeta = (int)floor((cotbeta - cotbetaLowEdge_) / cotbetaBinWidth_);
0478   iqbin = qbin > 2 ? 3 : qbin;
0479   if (icotalpha >= 0 && icotalpha < cotalphaBins_ && icotbeta >= 0 && icotbeta < cotbetaBins_) {
0480     qbinHist_[icotbeta][icotalpha]->Fill((double)iqbin);
0481     if (nxpix == 1)
0482       resSinglePixelXHist_[icotbeta][icotalpha]->Fill(dx / cmtomicron);
0483     else
0484       resMultiPixelXHist_[icotbeta][icotalpha][iqbin]->Fill(dx / cmtomicron);
0485     if (nypix == 1)
0486       resSinglePixelYHist_[icotbeta][icotalpha]->Fill(dy / cmtomicron);
0487     else
0488       resMultiPixelYHist_[icotbeta][icotalpha][iqbin]->Fill(dy / cmtomicron);
0489   }
0490
0491   return 0;
0492 }
0493
0494 //------------------------------------------------------------------------------
0495 //  Return the histogram generator for resolution in X.  A generator contains
0496 //  both the histogram and knows how to throw a random number off it.  It is
0497 //  called from FastSim (from PixelTemplateSmearerBase).
0498 //  If cotalpha or cotbeta are outside of the range, return the end of the range.
0499 //------------------------------------------------------------------------------
0500 const SimpleHistogramGenerator* PixelResolutionHistograms::getGeneratorX(double cotalpha,
0501                                                                          double cotbeta,
0502                                                                          int qbin,
0503                                                                          bool single) {
0504   int icotalpha, icotbeta, iqbin;
0505   icotalpha = (int)floor((cotalpha - cotalphaLowEdge_) / cotalphaBinWidth_);
0506   icotbeta = (int)floor((cotbeta - cotbetaLowEdge_) / cotbetaBinWidth_);
0507   iqbin = qbin > 2 ? 3 : qbin;  // if (qbin>2) then = 3, else return qbin
0508   //
0509   //if( icotalpha >= 0 && icotalpha < cotalphaBins_ && icotbeta >= 0 && icotbeta < cotbetaBins_ ) {
0510
0511   if (icotalpha < 0)
0512     icotalpha = 0;
0513   if (icotalpha >= cotalphaBins_)
0514     icotalpha = cotalphaBins_ - 1;
0515
0516   if (icotbeta < 0)
0517     icotbeta = 0;
0518   if (icotbeta >= cotbetaBins_)
0519     icotbeta = cotbetaBins_ - 1;
0520
0521   // At this point we are sure to return *some bin* from the 3D histogram
0522
0523   if (single)
0524     return resSinglePixelXGen_[icotbeta][icotalpha];
0525   else
0526     return resMultiPixelXGen_[icotbeta][icotalpha][iqbin];
0527
0528   // }
0529   //else
0530   //return nullptr;
0531 }
0532
0533 //------------------------------------------------------------------------------
0534 //  Return the histogram generator for resolution in Y.  A generator contains
0535 //  both the histogram and knows how to throw a random number off it.  It is
0536 //  called from FastSim (from PixelTemplateSmearerBase).
0537 //  If cotalpha or cotbeta are outside of the range, return the end of the range.
0538 //------------------------------------------------------------------------------
0539 const SimpleHistogramGenerator* PixelResolutionHistograms::getGeneratorY(double cotalpha,
0540                                                                          double cotbeta,
0541                                                                          int qbin,
0542                                                                          bool single) {
0543   int icotalpha, icotbeta, iqbin;
0544   icotalpha = (int)floor((cotalpha - cotalphaLowEdge_) / cotalphaBinWidth_);
0545   icotbeta = (int)floor((cotbeta - cotbetaLowEdge_) / cotbetaBinWidth_);
0546   iqbin = qbin > 2 ? 3 : qbin;  // if (qbin>2) then = 3, else return qbin
0547   //
0548   //if( icotalpha >= 0 && icotalpha < cotalphaBins_ && icotbeta >= 0 && icotbeta < cotbetaBins_ ) {
0549
0550   if (icotalpha < 0)
0551     icotalpha = 0;
0552   if (icotalpha >= cotalphaBins_)
0553     icotalpha = cotalphaBins_ - 1;
0554
0555   if (icotbeta < 0)
0556     icotbeta = 0;
0557   if (icotbeta >= cotbetaBins_)
0558     icotbeta = cotbetaBins_ - 1;
0559
0560   // At this point we are sure to return *some bin* from the 3D histogram
0561
0562   if (single)
0563     return resSinglePixelYGen_[icotbeta][icotalpha];
0564   else
0565     return resMultiPixelYGen_[icotbeta][icotalpha][iqbin];
0566
0567   //}
0568   //else
0569   //return nullptr;
0570 }