File indexing completed on 2024-09-07 04:36:19
0001
0002
0003
0004
0005
0006
0007
0008
0009 #ifdef SI_PIXEL_TEMPLATE_STANDALONE
0010
0011
0012
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_;
0024 };
0025 #define LOGDEBUG std::cout
0026 #define LOGERROR std::cout
0027 #define LOGINFO std::cout
0028
0029 #else
0030
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
0044 #include <cmath>
0045 #include <iostream>
0046 #include <string>
0047
0048
0049 #include <TFile.h>
0050 #include <TH1F.h>
0051 #include <TH2F.h>
0052
0053
0054 const float cmtomicron = 10000.0;
0055
0056
0057
0058
0059
0060
0061
0062 PixelResolutionHistograms::PixelResolutionHistograms(std::string filename,
0063 std::string rootdir,
0064 std::string descTitle,
0065 unsigned int detType,
0066 double cotbetaBinWidth,
0067 double cotbetaLowEdge,
0068 int cotbetaBins,
0069 double cotalphaBinWidth,
0070 double cotalphaLowEdge,
0071 int cotalphaBins)
0072
0073
0074 : weOwnHistograms_(true),
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_(),
0087 resMultiPixelYHist_(),
0088 resSinglePixelYHist_(),
0089 qbinHist_(),
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
0099
0100
0101
0102
0103
0104
0105
0106
0107
0108
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
0119 binningHisto_->SetBinContent(0, 0, detType_);
0120
0121
0122 Char_t histo[200];
0123 Char_t title[200];
0124
0125
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);
0132
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
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);
0163
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
0198
0199
0200
0201
0202
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),
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_(),
0223 resMultiPixelYHist_(),
0224 resSinglePixelYHist_(),
0225 qbinHist_(),
0226 file_(nullptr),
0227 status_(0),
0228 resMultiPixelXGen_(),
0229 resSinglePixelXGen_(),
0230 resMultiPixelYGen_(),
0231 resSinglePixelYGen_(),
0232 qbinGen_() {
0233 Char_t histo[200];
0234 Char_t title[200];
0235 TH1F* tmphist = nullptr;
0236
0237
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;
0243 }
0244
0245
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;
0251 }
0252
0253 if (detType == -1) {
0254
0255 detType_ = binningHisto_->GetBinContent(0, 0);
0256 } else {
0257 detType_ = detType;
0258 }
0259
0260
0261 cotbetaAxis_ = binningHisto_->GetXaxis();
0262 cotbetaBinWidth_ = binningHisto_->GetXaxis()->GetBinWidth(1);
0263 cotbetaLowEdge_ = binningHisto_->GetXaxis()->GetXmin();
0264 cotbetaBins_ = binningHisto_->GetXaxis()->GetNbins();
0265 cotalphaAxis_ = binningHisto_->GetYaxis();
0266 cotalphaBinWidth_ = binningHisto_->GetYaxis()->GetBinWidth(1);
0267 cotalphaLowEdge_ = binningHisto_->GetYaxis()->GetXmin();
0268 cotalphaBins_ = binningHisto_->GetYaxis()->GetNbins();
0269
0270 if (!ignore_multi) {
0271
0272
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);
0279
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 }
0331
0332
0333
0334
0335 for (int ii = 0; ii < cotbetaBins_; ii++) {
0336 for (int jj = 0; jj < cotalphaBins_; jj++) {
0337
0338
0339 sprintf(histo, "hx%d0%02d%d", detType_, ii + 1, jj + 1);
0340
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
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
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
0428
0429
0430
0431 PixelResolutionHistograms::~PixelResolutionHistograms() {
0432
0433
0434
0435 if (!weOwnHistograms_) {
0436
0437 file_->Close();
0438
0439
0440 } else {
0441
0442 LOGINFO << "PixelResHistoStore: Writing the histograms to the output file. "
0443 << std::endl;
0444 file_->Write();
0445 file_->Close();
0446
0447
0448
0449 }
0450
0451
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
0471
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
0496
0497
0498
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;
0508
0509
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
0522
0523 if (single)
0524 return resSinglePixelXGen_[icotbeta][icotalpha];
0525 else
0526 return resMultiPixelXGen_[icotbeta][icotalpha][iqbin];
0527
0528
0529
0530
0531 }
0532
0533
0534
0535
0536
0537
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;
0547
0548
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
0561
0562 if (single)
0563 return resSinglePixelYGen_[icotbeta][icotalpha];
0564 else
0565 return resMultiPixelYGen_[icotbeta][icotalpha][iqbin];
0566
0567
0568
0569
0570 }