File indexing completed on 2024-04-06 11:59:03
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020 #include <memory>
0021 #include <string>
0022 #include <iostream>
0023 #include <fstream>
0024 #include <sstream>
0025 #include <vector>
0026 #include <map>
0027
0028
0029 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0030 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0031 #include "FWCore/Common/interface/TriggerNames.h"
0032 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0033 #include "FWCore/Framework/interface/Event.h"
0034 #include "FWCore/Framework/interface/EventSetup.h"
0035 #include "FWCore/Framework/interface/Frameworkfwd.h"
0036 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0037 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0038 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0039 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0040 #include "FWCore/ServiceRegistry/interface/Service.h"
0041 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0042 #include "DataFormats/Common/interface/TriggerResults.h"
0043 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0044 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0045 #include "DataFormats/HcalDigi/interface/HcalCalibrationEventTypes.h"
0046 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0047 #include "DataFormats/HcalDigi/interface/HBHEDataFrame.h"
0048 #include "DataFormats/HcalDigi/interface/HFDataFrame.h"
0049 #include "DataFormats/HcalDigi/interface/HODataFrame.h"
0050 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0051 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
0052 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
0053 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0054 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0055
0056 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0057
0058 #include "TH1D.h"
0059 #include "TH2D.h"
0060 #include "TMath.h"
0061 #include "TProfile.h"
0062 #include "TTree.h"
0063
0064
0065
0066
0067 class RecAnalyzerMinbias : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0068 public:
0069 explicit RecAnalyzerMinbias(const edm::ParameterSet&);
0070 ~RecAnalyzerMinbias() override = default;
0071
0072 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0073
0074 private:
0075 void analyze(edm::Event const&, edm::EventSetup const&) override;
0076 void beginJob() override;
0077 void endJob() override;
0078 void beginRun(edm::Run const&, edm::EventSetup const&) override;
0079 void endRun(edm::Run const&, edm::EventSetup const&) override {}
0080
0081 private:
0082 void analyzeHcal(const HBHERecHitCollection&, const HFRecHitCollection&, int, bool, double);
0083
0084
0085 const bool runNZS_, Noise_, ignoreL1_;
0086 const bool fillHist_, extraHist_;
0087 const double eLowHB_, eHighHB_, eLowHE_, eHighHE_;
0088 const double eLowHF_, eHighHF_, eMin_;
0089 const int runMin_, runMax_;
0090 const std::vector<int> trigbit_;
0091 const std::string cfile_;
0092 bool theRecalib_, init_;
0093 std::map<DetId, double> corrFactor_;
0094 std::vector<unsigned int> hcalID_;
0095 TTree *myTree_, *myTree1_;
0096 TH1D* h_[4];
0097 TH2D *hbhe_, *hb_, *he_, *hf_;
0098 TH1D *h_AmplitudeHBtest_, *h_AmplitudeHEtest_;
0099 TH1D* h_AmplitudeHFtest_;
0100 TH1D *h_AmplitudeHB_, *h_AmplitudeHE_, *h_AmplitudeHF_;
0101 TProfile *hbherun_, *hbrun_, *herun_, *hfrun_;
0102 std::vector<TH1D*> histo_;
0103 std::map<HcalDetId, TH1D*> histHC_;
0104 double rnnum_;
0105 struct myInfo {
0106 double theMB0, theMB1, theMB2, theMB3, theMB4, runcheck;
0107 myInfo() { theMB0 = theMB1 = theMB2 = theMB3 = theMB4 = runcheck = 0; }
0108 };
0109
0110 double rnnumber;
0111 int mysubd, depth, iphi, ieta, cells, trigbit;
0112 float mom0_MB, mom1_MB, mom2_MB, mom3_MB, mom4_MB;
0113 int HBHEsize, HFsize;
0114 std::map<std::pair<int, HcalDetId>, myInfo> myMap_;
0115 const edm::EDGetTokenT<HBHERecHitCollection> tok_hbherecoMB_;
0116 const edm::EDGetTokenT<HFRecHitCollection> tok_hfrecoMB_;
0117 const edm::EDGetTokenT<L1GlobalTriggerObjectMapRecord> tok_hltL1GtMap_;
0118 const edm::EDGetTokenT<GenEventInfoProduct> tok_ew_;
0119 const edm::EDGetTokenT<HBHEDigiCollection> tok_hbhedigi_;
0120 const edm::EDGetTokenT<QIE11DigiCollection> tok_qie11digi_;
0121 const edm::EDGetTokenT<HODigiCollection> tok_hodigi_;
0122 const edm::EDGetTokenT<HFDigiCollection> tok_hfdigi_;
0123 const edm::EDGetTokenT<QIE10DigiCollection> tok_qie10digi_;
0124 const edm::EDGetTokenT<L1GlobalTriggerReadoutRecord> tok_gtRec_;
0125 const edm::ESGetToken<HcalTopology, HcalRecNumberingRecord> tok_htopo_;
0126 };
0127
0128
0129 RecAnalyzerMinbias::RecAnalyzerMinbias(const edm::ParameterSet& iConfig)
0130 : runNZS_(iConfig.getParameter<bool>("runNZS")),
0131 Noise_(iConfig.getParameter<bool>("noise")),
0132 ignoreL1_(iConfig.getUntrackedParameter<bool>("ignoreL1", false)),
0133 fillHist_(iConfig.getUntrackedParameter<bool>("fillHisto", false)),
0134 extraHist_(iConfig.getUntrackedParameter<bool>("extraHisto", false)),
0135 eLowHB_(iConfig.getParameter<double>("eLowHB")),
0136 eHighHB_(iConfig.getParameter<double>("eHighHB")),
0137 eLowHE_(iConfig.getParameter<double>("eLowHE")),
0138 eHighHE_(iConfig.getParameter<double>("eHighHE")),
0139 eLowHF_(iConfig.getParameter<double>("eLowHF")),
0140 eHighHF_(iConfig.getParameter<double>("eHighHF")),
0141 eMin_(iConfig.getUntrackedParameter<double>("eMin", 2.0)),
0142 runMin_(iConfig.getUntrackedParameter<int>("RunMin", 308327)),
0143 runMax_(iConfig.getUntrackedParameter<int>("RunMax", 315250)),
0144 trigbit_(iConfig.getUntrackedParameter<std::vector<int>>("triggerBits")),
0145 cfile_(iConfig.getUntrackedParameter<std::string>("corrFile")),
0146 init_(false),
0147 tok_hbherecoMB_(consumes<HBHERecHitCollection>(iConfig.getParameter<edm::InputTag>("hbheInputMB"))),
0148 tok_hfrecoMB_(consumes<HFRecHitCollection>(iConfig.getParameter<edm::InputTag>("hfInputMB"))),
0149 tok_hltL1GtMap_(consumes<L1GlobalTriggerObjectMapRecord>(edm::InputTag("hltL1GtObjectMap"))),
0150 tok_ew_(consumes<GenEventInfoProduct>(edm::InputTag("generator"))),
0151 tok_hbhedigi_(consumes<HBHEDigiCollection>(iConfig.getParameter<edm::InputTag>("hcalDigiCollectionTag"))),
0152 tok_qie11digi_(consumes<QIE11DigiCollection>(iConfig.getParameter<edm::InputTag>("hcalDigiCollectionTag"))),
0153 tok_hodigi_(consumes<HODigiCollection>(iConfig.getParameter<edm::InputTag>("hcalDigiCollectionTag"))),
0154 tok_hfdigi_(consumes<HFDigiCollection>(iConfig.getParameter<edm::InputTag>("hcalDigiCollectionTag"))),
0155 tok_qie10digi_(consumes<QIE10DigiCollection>(iConfig.getParameter<edm::InputTag>("hcalDigiCollectionTag"))),
0156 tok_gtRec_(consumes<L1GlobalTriggerReadoutRecord>(edm::InputTag("gtDigisAlCaMB"))),
0157 tok_htopo_(esConsumes<HcalTopology, HcalRecNumberingRecord, edm::Transition::BeginRun>()) {
0158 usesResource("TFileService");
0159
0160
0161 std::vector<int> ieta = iConfig.getUntrackedParameter<std::vector<int>>("hcalIeta");
0162 std::vector<int> iphi = iConfig.getUntrackedParameter<std::vector<int>>("hcalIphi");
0163 std::vector<int> depth = iConfig.getUntrackedParameter<std::vector<int>>("hcalDepth");
0164
0165
0166 std::ifstream infile(cfile_.c_str());
0167 if (!infile.is_open()) {
0168 theRecalib_ = false;
0169 edm::LogWarning("RecAnalyzerMinbias") << "Cannot open '" << cfile_ << "' for the correction file";
0170 } else {
0171 unsigned int ndets(0), nrec(0);
0172 while (true) {
0173 unsigned int id;
0174 double cfac;
0175 infile >> id >> cfac;
0176 if (!infile.good())
0177 break;
0178 HcalDetId detId(id);
0179 nrec++;
0180 std::map<DetId, double>::iterator itr = corrFactor_.find(detId);
0181 if (itr == corrFactor_.end()) {
0182 corrFactor_[detId] = cfac;
0183 ndets++;
0184 }
0185 }
0186 infile.close();
0187 edm::LogVerbatim("RecAnalyzerMinbias") << "Reads " << nrec << " correction factors for " << ndets << " detIds";
0188 theRecalib_ = (ndets > 0);
0189 }
0190
0191 edm::LogVerbatim("RecAnalyzerMinbias") << " Flags (ReCalib): " << theRecalib_ << " (IgnoreL1): " << ignoreL1_
0192 << " (NZS) " << runNZS_ << " and with " << ieta.size()
0193 << " detId for full histogram";
0194 edm::LogVerbatim("RecAnalyzerMinbias") << "Thresholds for HB " << eLowHB_ << ":" << eHighHB_ << " for HE " << eLowHE_
0195 << ":" << eHighHE_ << " for HF " << eLowHF_ << ":" << eHighHF_;
0196 for (unsigned int k = 0; k < ieta.size(); ++k) {
0197 HcalSubdetector subd = ((std::abs(ieta[k]) > 29) ? HcalForward
0198 : (std::abs(ieta[k]) > 16) ? HcalEndcap
0199 : ((std::abs(ieta[k]) == 16) && (depth[k] == 3)) ? HcalEndcap
0200 : (depth[k] == 4) ? HcalOuter
0201 : HcalBarrel);
0202 unsigned int id = (HcalDetId(subd, ieta[k], iphi[k], depth[k])).rawId();
0203 hcalID_.push_back(id);
0204 edm::LogVerbatim("RecAnalyzerMinbias") << "DetId[" << k << "] " << HcalDetId(id);
0205 }
0206 edm::LogVerbatim("RecAnalyzerMinbias") << "Select on " << trigbit_.size() << " L1 Trigger selection";
0207 for (unsigned int k = 0; k < trigbit_.size(); ++k)
0208 edm::LogVerbatim("RecAnalyzerMinbias") << "Bit[" << k << "] " << trigbit_[k];
0209 }
0210
0211 void RecAnalyzerMinbias::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0212 std::vector<int> iarray;
0213 edm::ParameterSetDescription desc;
0214 desc.add<bool>("runNZS", true);
0215 desc.add<bool>("noise", false);
0216 desc.add<double>("eLowHB", 4);
0217 desc.add<double>("eHighHB", 100);
0218 desc.add<double>("eLowHE", 4);
0219 desc.add<double>("eHighHE", 150);
0220 desc.add<double>("eLowHF", 10);
0221 desc.add<double>("eHighHF", 150);
0222
0223 desc.addUntracked<double>("eMin", 2.0);
0224
0225 desc.addUntracked<int>("runMin", 308327);
0226 desc.addUntracked<int>("runMax", 308347);
0227 desc.addUntracked<std::vector<int>>("triggerBits", iarray);
0228 desc.addUntracked<bool>("ignoreL1", false);
0229 desc.addUntracked<std::string>("corrFile", "CorFactor.txt");
0230 desc.addUntracked<bool>("fillHisto", false);
0231 desc.addUntracked<bool>("extraHisto", false);
0232 desc.addUntracked<std::vector<int>>("hcalIeta", iarray);
0233 desc.addUntracked<std::vector<int>>("hcalIphi", iarray);
0234 desc.addUntracked<std::vector<int>>("hcalDepth", iarray);
0235 desc.add<edm::InputTag>("hbheInputMB", edm::InputTag("hbherecoMB"));
0236 desc.add<edm::InputTag>("hfInputMB", edm::InputTag("hfrecoMB"));
0237 desc.add<edm::InputTag>("gtDigisAlCaMB", edm::InputTag("gtDigisAlCaMB"));
0238 desc.add<edm::InputTag>("hcalDigiCollectionTag", edm::InputTag("hcalDigis"));
0239 descriptions.add("recAnalyzerMinbias", desc);
0240 }
0241
0242 void RecAnalyzerMinbias::beginJob() {
0243 edm::Service<TFileService> fs;
0244 std::string hc[5] = {"Empty", "HB", "HE", "HO", "HF"};
0245 char name[700], title[700];
0246 hbhe_ = fs->make<TH2D>("hbhe", "Noise in HB/HE", 61, -30.5, 30.5, 72, 0.5, 72.5);
0247 hb_ = fs->make<TH2D>("hb", "Noise in HB", 61, -16.5, 16.5, 72, 0.5, 72.5);
0248 he_ = fs->make<TH2D>("he", "Noise in HE", 61, -30.5, 30.5, 72, 0.5, 72.5);
0249 hf_ = fs->make<TH2D>("hf", "Noise in HF", 82, -41.5, 41.5, 72, 0.5, 72.5);
0250 int nbin = (runMax_ - runMin_ + 1);
0251 sprintf(title, "Fraction of channels in HB/HE with E > %4.1f GeV vs Run number", eMin_);
0252 hbherun_ = fs->make<TProfile>("hbherun", title, nbin, runMin_ - 0.5, runMax_ + 0.5, 0.0, 1.0);
0253 sprintf(title, "Fraction of channels in HB with E > %4.1f GeV vs Run number", eMin_);
0254 hbrun_ = fs->make<TProfile>("hbrun", title, nbin, runMin_ - 0.5, runMax_ + 0.5, 0.0, 1.0);
0255 sprintf(title, "Fraction of channels in HE with E > %4.1f GeV vs Run number", eMin_);
0256 herun_ = fs->make<TProfile>("herun", title, nbin, runMin_ - 0.5, runMax_ + 0.5, 0.0, 1.0);
0257 sprintf(title, "Fraction of channels in HF with E > %4.1f GeV vs Run number", eMin_);
0258 hfrun_ = fs->make<TProfile>("hfrun", title, nbin, runMin_ - 0.5, runMax_ + 0.5, 0.0, 1.0);
0259 for (int idet = 1; idet <= 4; idet++) {
0260 sprintf(name, "%s", hc[idet].c_str());
0261 sprintf(title, "Noise distribution for %s", hc[idet].c_str());
0262 h_[idet - 1] = fs->make<TH1D>(name, title, 48, -6., 6.);
0263 }
0264
0265 for (const auto& hcalid : hcalID_) {
0266 HcalDetId id = HcalDetId(hcalid);
0267 int subdet = id.subdetId();
0268 sprintf(name, "%s%d_%d_%d", hc[subdet].c_str(), id.ieta(), id.iphi(), id.depth());
0269 sprintf(title,
0270 "Energy Distribution for %s ieta %d iphi %d depth %d",
0271 hc[subdet].c_str(),
0272 id.ieta(),
0273 id.iphi(),
0274 id.depth());
0275 double xmin = (subdet == 4) ? -10 : -1;
0276 double xmax = (subdet == 4) ? 90 : 9;
0277 TH1D* hh = fs->make<TH1D>(name, title, 50, xmin, xmax);
0278 histo_.push_back(hh);
0279 };
0280
0281 if (extraHist_) {
0282 h_AmplitudeHBtest_ = fs->make<TH1D>("h_AmplitudeHBtest", "", 5000, 0., 5000.);
0283 h_AmplitudeHEtest_ = fs->make<TH1D>("h_AmplitudeHEtest", "", 3000, 0., 3000.);
0284 h_AmplitudeHFtest_ = fs->make<TH1D>("h_AmplitudeHFtest", "", 10000, 0., 10000.);
0285 h_AmplitudeHB_ = fs->make<TH1D>("h_AmplitudeHB", "", 100000, 0., 100000.);
0286 h_AmplitudeHE_ = fs->make<TH1D>("h_AmplitudeHE", "", 300000, 0., 300000.);
0287 h_AmplitudeHF_ = fs->make<TH1D>("h_AmplitudeHF", "", 100000, 0., 1000000.);
0288 }
0289
0290 if (!fillHist_) {
0291 myTree_ = fs->make<TTree>("RecJet", "RecJet Tree");
0292 myTree_->Branch("cells", &cells, "cells/I");
0293 myTree_->Branch("mysubd", &mysubd, "mysubd/I");
0294 myTree_->Branch("depth", &depth, "depth/I");
0295 myTree_->Branch("ieta", &ieta, "ieta/I");
0296 myTree_->Branch("iphi", &iphi, "iphi/I");
0297 myTree_->Branch("mom0_MB", &mom0_MB, "mom0_MB/F");
0298 myTree_->Branch("mom1_MB", &mom1_MB, "mom1_MB/F");
0299 myTree_->Branch("mom2_MB", &mom2_MB, "mom2_MB/F");
0300 myTree_->Branch("mom3_MB", &mom3_MB, "mom3_MB/F");
0301 myTree_->Branch("mom4_MB", &mom4_MB, "mom4_MB/F");
0302 myTree_->Branch("trigbit", &trigbit, "trigbit/I");
0303 myTree_->Branch("rnnumber", &rnnumber, "rnnumber/D");
0304 }
0305 myTree1_ = fs->make<TTree>("RecJet1", "RecJet1 Tree");
0306 myTree1_->Branch("rnnum_", &rnnum_, "rnnum_/D");
0307 myTree1_->Branch("HBHEsize", &HBHEsize, "HBHEsize/I");
0308 myTree1_->Branch("HFsize", &HFsize, "HFsize/I");
0309
0310 myMap_.clear();
0311 }
0312
0313
0314
0315 void RecAnalyzerMinbias::endJob() {
0316 if (!fillHist_) {
0317 cells = 0;
0318 for (const auto& itr : myMap_) {
0319 edm::LogVerbatim("RecAnalyzerMinbias") << "Fired trigger bit number " << itr.first.first;
0320 myInfo info = itr.second;
0321 if (info.theMB0 > 0) {
0322 mom0_MB = info.theMB0;
0323 mom1_MB = info.theMB1;
0324 mom2_MB = info.theMB2;
0325 mom3_MB = info.theMB3;
0326 mom4_MB = info.theMB4;
0327 rnnumber = info.runcheck;
0328 trigbit = itr.first.first;
0329 mysubd = itr.first.second.subdet();
0330 depth = itr.first.second.depth();
0331 iphi = itr.first.second.iphi();
0332 ieta = itr.first.second.ieta();
0333 edm::LogVerbatim("RecAnalyzerMinbias")
0334 << " Result= " << trigbit << " " << mysubd << " " << ieta << " " << iphi << " mom0 " << mom0_MB
0335 << " mom1 " << mom1_MB << " mom2 " << mom2_MB << " mom3 " << mom3_MB << " mom4 " << mom4_MB;
0336 myTree_->Fill();
0337 cells++;
0338 }
0339 }
0340 edm::LogVerbatim("RecAnalyzerMinbias") << "cells " << cells;
0341 }
0342 #ifdef EDM_ML_DEBUG
0343 edm::LogVerbatim("RecAnalyzerMinbias") << "Exiting from RecAnalyzerMinbias::endjob";
0344 #endif
0345 }
0346
0347 void RecAnalyzerMinbias::beginRun(edm::Run const&, edm::EventSetup const& iS) {
0348 if (!init_) {
0349 init_ = true;
0350 if (fillHist_) {
0351 edm::Service<TFileService> fs;
0352 const HcalTopology* hcaltopology = &iS.getData(tok_htopo_);
0353
0354 char name[700], title[700];
0355
0356 int maxDepthHB = hcaltopology->maxDepthHB();
0357 int nbinHB = (Noise_) ? 18 : int(2000 * eHighHB_);
0358 double x_min = (Noise_) ? -3. : 0.;
0359 double x_max = (Noise_) ? 3. : 2. * eHighHB_;
0360 for (int eta = -50; eta < 50; eta++) {
0361 for (int phi = 0; phi < 100; phi++) {
0362 for (int depth = 1; depth <= maxDepthHB; depth++) {
0363 HcalDetId cell(HcalBarrel, eta, phi, depth);
0364 if (hcaltopology->valid(cell)) {
0365 sprintf(name, "HBeta%dphi%ddep%d", eta, phi, depth);
0366 sprintf(title, "HB #eta %d #phi %d depth %d", eta, phi, depth);
0367 TH1D* h = fs->make<TH1D>(name, title, nbinHB, x_min, x_max);
0368 histHC_[cell] = h;
0369 }
0370 }
0371 }
0372 }
0373
0374 int maxDepthHE = hcaltopology->maxDepthHE();
0375 int nbinHE = (Noise_) ? 18 : int(2000 * eHighHE_);
0376 x_min = (Noise_) ? -3. : 0.;
0377 x_max = (Noise_) ? 3. : 2. * eHighHE_;
0378 for (int eta = -50; eta < 50; eta++) {
0379 for (int phi = 0; phi < 100; phi++) {
0380 for (int depth = 1; depth <= maxDepthHE; depth++) {
0381 HcalDetId cell(HcalEndcap, eta, phi, depth);
0382 if (hcaltopology->valid(cell)) {
0383 sprintf(name, "HEeta%dphi%ddep%d", eta, phi, depth);
0384 sprintf(title, "HE #eta %d #phi %d depth %d", eta, phi, depth);
0385 TH1D* h = fs->make<TH1D>(name, title, nbinHE, x_min, x_max);
0386 histHC_[cell] = h;
0387 }
0388 }
0389 }
0390 }
0391
0392 int maxDepthHF = 4;
0393 int nbinHF = (Noise_) ? 200 : int(2000 * eHighHF_);
0394 x_min = (Noise_) ? -10. : 0.;
0395 x_max = (Noise_) ? 10. : 2. * eHighHF_;
0396 for (int eta = -50; eta < 50; eta++) {
0397 for (int phi = 0; phi < 100; phi++) {
0398 for (int depth = 1; depth <= maxDepthHF; depth++) {
0399 HcalDetId cell(HcalForward, eta, phi, depth);
0400 if (hcaltopology->valid(cell)) {
0401 sprintf(name, "HFeta%dphi%ddep%d", eta, phi, depth);
0402 sprintf(title, "Energy (HF #eta %d #phi %d depth %d)", eta, phi, depth);
0403 TH1D* h = fs->make<TH1D>(name, title, nbinHF, x_min, x_max);
0404 histHC_[cell] = h;
0405 }
0406 }
0407 }
0408 }
0409 }
0410 }
0411 }
0412
0413
0414
0415
0416
0417
0418 void RecAnalyzerMinbias::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0419 rnnum_ = (double)iEvent.run();
0420
0421 if (extraHist_) {
0422 double amplitudefullHB(0), amplitudefullHE(0), amplitudefullHF(0);
0423 const edm::Handle<HBHEDigiCollection>& hbhedigi = iEvent.getHandle(tok_hbhedigi_);
0424 if (hbhedigi.isValid()) {
0425 for (auto const& digi : *(hbhedigi.product())) {
0426 int nTS = digi.size();
0427 double amplitudefullTSs = 0.;
0428 if (digi.id().subdet() == HcalBarrel) {
0429 if (nTS <= 10) {
0430 for (int i = 0; i < nTS; i++)
0431 amplitudefullTSs += digi.sample(i).adc();
0432 h_AmplitudeHBtest_->Fill(amplitudefullTSs);
0433 amplitudefullHB += amplitudefullTSs;
0434 }
0435 }
0436 if (digi.id().subdet() == HcalEndcap) {
0437 if (nTS <= 10) {
0438 for (int i = 0; i < nTS; i++)
0439 amplitudefullTSs += digi.sample(i).adc();
0440 h_AmplitudeHEtest_->Fill(amplitudefullTSs);
0441 amplitudefullHE += amplitudefullTSs;
0442 }
0443 }
0444 }
0445 }
0446
0447 const edm::Handle<QIE11DigiCollection>& qie11digi = iEvent.getHandle(tok_qie11digi_);
0448 if (qie11digi.isValid()) {
0449 for (QIE11DataFrame const digi : *(qie11digi.product())) {
0450 double amplitudefullTSs = 0.;
0451 if (HcalDetId(digi.id()).subdet() == HcalBarrel) {
0452 for (int i = 0; i < digi.samples(); i++)
0453 amplitudefullTSs += digi[i].adc();
0454 h_AmplitudeHBtest_->Fill(amplitudefullTSs);
0455 amplitudefullHB += amplitudefullTSs;
0456 }
0457 if (HcalDetId(digi.id()).subdet() == HcalEndcap) {
0458 for (int i = 0; i < digi.samples(); i++)
0459 amplitudefullTSs += digi[i].adc();
0460 h_AmplitudeHEtest_->Fill(amplitudefullTSs);
0461 amplitudefullHE += amplitudefullTSs;
0462 }
0463 }
0464 }
0465
0466 const edm::Handle<HFDigiCollection>& hfdigi = iEvent.getHandle(tok_hfdigi_);
0467 if (hfdigi.isValid()) {
0468 for (auto const& digi : *(hfdigi.product())) {
0469 int nTS = digi.size();
0470 double amplitudefullTSs = 0.;
0471 if (digi.id().subdet() == HcalForward) {
0472 if (nTS <= 10) {
0473 for (int i = 0; i < nTS; i++)
0474 amplitudefullTSs += digi.sample(i).adc();
0475 h_AmplitudeHFtest_->Fill(amplitudefullTSs);
0476 amplitudefullHF += amplitudefullTSs;
0477 }
0478 }
0479 }
0480 }
0481
0482 const edm::Handle<QIE10DigiCollection>& qie10digi = iEvent.getHandle(tok_qie10digi_);
0483 if (qie10digi.isValid()) {
0484 for (QIE10DataFrame const digi : *(qie10digi.product())) {
0485 double amplitudefullTSs = 0.;
0486 if (HcalDetId(digi.id()).subdet() == HcalForward) {
0487 for (int i = 0; i < digi.samples(); i++)
0488 amplitudefullTSs += digi[i].adc();
0489 h_AmplitudeHFtest_->Fill(amplitudefullTSs);
0490 amplitudefullHF += amplitudefullTSs;
0491 }
0492 }
0493 }
0494
0495 h_AmplitudeHB_->Fill(amplitudefullHB);
0496 h_AmplitudeHE_->Fill(amplitudefullHE);
0497 h_AmplitudeHF_->Fill(amplitudefullHF);
0498 }
0499
0500 const edm::Handle<HBHERecHitCollection>& hbheMB = iEvent.getHandle(tok_hbherecoMB_);
0501 if (!hbheMB.isValid()) {
0502 edm::LogWarning("RecAnalyzerMinbias") << "HcalCalibAlgos: Error! can't get hbhe product!";
0503 return;
0504 }
0505 const HBHERecHitCollection HithbheMB = *(hbheMB.product());
0506 HBHEsize = HithbheMB.size();
0507 edm::LogVerbatim("RecAnalyzerMinbias") << "HBHE MB size of collection " << HithbheMB.size();
0508 if (HithbheMB.size() < 5100 && runNZS_) {
0509 edm::LogWarning("RecAnalyzerMinbias") << "HBHE problem " << rnnum_ << " size " << HBHEsize;
0510 }
0511
0512 const edm::Handle<HFRecHitCollection> hfMB = iEvent.getHandle(tok_hfrecoMB_);
0513 if (!hfMB.isValid()) {
0514 edm::LogWarning("RecAnalyzerMinbias") << "HcalCalibAlgos: Error! can't get hf product!";
0515 return;
0516 }
0517 const HFRecHitCollection HithfMB = *(hfMB.product());
0518 edm::LogVerbatim("RecAnalyzerMinbias") << "HF MB size of collection " << HithfMB.size();
0519 HFsize = HithfMB.size();
0520 if (HithfMB.size() < 1700 && runNZS_) {
0521 edm::LogWarning("RecAnalyzerMinbias") << "HF problem " << rnnum_ << " size " << HFsize;
0522 }
0523
0524 bool select(false);
0525 if (!trigbit_.empty()) {
0526 const edm::Handle<L1GlobalTriggerObjectMapRecord>& gtObjectMapRecord = iEvent.getHandle(tok_hltL1GtMap_);
0527 if (gtObjectMapRecord.isValid()) {
0528 const std::vector<L1GlobalTriggerObjectMap>& objMapVec = gtObjectMapRecord->gtObjectMap();
0529 for (std::vector<L1GlobalTriggerObjectMap>::const_iterator itMap = objMapVec.begin(); itMap != objMapVec.end();
0530 ++itMap) {
0531 bool resultGt = (*itMap).algoGtlResult();
0532 if (resultGt) {
0533 int algoBit = (*itMap).algoBitNumber();
0534 if (std::find(trigbit_.begin(), trigbit_.end(), algoBit) != trigbit_.end()) {
0535 select = true;
0536 break;
0537 }
0538 }
0539 }
0540 }
0541 }
0542
0543 if (!trigbit_.empty() || select)
0544 myTree1_->Fill();
0545
0546
0547 double eventWeight = 1.0;
0548 const edm::Handle<GenEventInfoProduct>& genEventInfo = iEvent.getHandle(tok_ew_);
0549 if (genEventInfo.isValid())
0550 eventWeight = genEventInfo->weight();
0551 #ifdef EDM_ML_DEBUG
0552 edm::LogVerbatim("RecAnalyzerMinbias") << "Test HB " << HBHEsize << " HF " << HFsize << " Trigger " << trigbit_.size()
0553 << ":" << select << ":" << ignoreL1_ << " Wt " << eventWeight;
0554 #endif
0555 if (ignoreL1_ || (!trigbit_.empty() && select)) {
0556 analyzeHcal(HithbheMB, HithfMB, 1, true, eventWeight);
0557 } else if ((!ignoreL1_) && (trigbit_.empty())) {
0558 const edm::Handle<L1GlobalTriggerObjectMapRecord>& gtObjectMapRecord = iEvent.getHandle(tok_hltL1GtMap_);
0559 if (gtObjectMapRecord.isValid()) {
0560 const std::vector<L1GlobalTriggerObjectMap>& objMapVec = gtObjectMapRecord->gtObjectMap();
0561 bool ok(false);
0562 for (std::vector<L1GlobalTriggerObjectMap>::const_iterator itMap = objMapVec.begin(); itMap != objMapVec.end();
0563 ++itMap) {
0564 bool resultGt = (*itMap).algoGtlResult();
0565 if (resultGt) {
0566 int algoBit = (*itMap).algoBitNumber();
0567 analyzeHcal(HithbheMB, HithfMB, algoBit, (!ok), eventWeight);
0568 ok = true;
0569 }
0570 }
0571 if (!ok) {
0572 edm::LogVerbatim("RecAnalyzerMinbias") << "No passed L1 Trigger found";
0573 }
0574 }
0575 }
0576 }
0577
0578 void RecAnalyzerMinbias::analyzeHcal(
0579 const HBHERecHitCollection& HithbheMB, const HFRecHitCollection& HithfMB, int algoBit, bool fill, double weight) {
0580
0581 int count(0), countHB(0), countHE(0), count2(0), count2HB(0), count2HE(0);
0582 for (HBHERecHitCollection::const_iterator hbheItr = HithbheMB.begin(); hbheItr != HithbheMB.end(); hbheItr++) {
0583
0584 DetId mydetid = hbheItr->id().rawId();
0585 double icalconst(1.);
0586 if (theRecalib_) {
0587 std::map<DetId, double>::iterator itr = corrFactor_.find(mydetid);
0588 if (itr != corrFactor_.end())
0589 icalconst = itr->second;
0590 }
0591 HBHERecHit aHit(hbheItr->id(), hbheItr->energy() * icalconst, hbheItr->time());
0592 double energyhit = aHit.energy();
0593 DetId id = (*hbheItr).detid();
0594 HcalDetId hid = HcalDetId(id);
0595 double eLow = (hid.subdet() == HcalEndcap) ? eLowHE_ : eLowHB_;
0596 double eHigh = (hid.subdet() == HcalEndcap) ? eHighHE_ : eHighHB_;
0597 ++count;
0598 if (id.subdetId() == HcalBarrel)
0599 ++countHB;
0600 else
0601 ++countHE;
0602 if (fill) {
0603 for (unsigned int i = 0; i < hcalID_.size(); i++) {
0604 if (hcalID_[i] == id.rawId()) {
0605 histo_[i]->Fill(energyhit);
0606 break;
0607 }
0608 }
0609 if (fillHist_) {
0610 std::map<HcalDetId, TH1D*>::iterator itr1 = histHC_.find(hid);
0611 if (itr1 != histHC_.end())
0612 itr1->second->Fill(energyhit);
0613 }
0614 h_[hid.subdet() - 1]->Fill(energyhit);
0615 if (energyhit > eMin_) {
0616 hbhe_->Fill(hid.ieta(), hid.iphi());
0617 ++count2;
0618 if (id.subdetId() == HcalBarrel) {
0619 ++count2HB;
0620 hb_->Fill(hid.ieta(), hid.iphi());
0621 } else {
0622 ++count2HE;
0623 he_->Fill(hid.ieta(), hid.iphi());
0624 }
0625 }
0626 }
0627 if (!fillHist_) {
0628 if (Noise_ || runNZS_ || (energyhit >= eLow && energyhit <= eHigh)) {
0629 std::map<std::pair<int, HcalDetId>, myInfo>::iterator itr1 =
0630 myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0631 if (itr1 == myMap_.end()) {
0632 myInfo info;
0633 myMap_[std::pair<int, HcalDetId>(algoBit, hid)] = info;
0634 itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0635 }
0636 itr1->second.theMB0 += weight;
0637 itr1->second.theMB1 += (weight * energyhit);
0638 itr1->second.theMB2 += (weight * energyhit * energyhit);
0639 itr1->second.theMB3 += (weight * energyhit * energyhit * energyhit);
0640 itr1->second.theMB4 += (weight * energyhit * energyhit * energyhit * energyhit);
0641 itr1->second.runcheck = rnnum_;
0642 }
0643 }
0644 }
0645 if (fill) {
0646 if (count > 0)
0647 hbherun_->Fill(rnnum_, (double)(count2) / count);
0648 if (countHB > 0)
0649 hbrun_->Fill(rnnum_, (double)(count2HB) / countHB);
0650 if (countHE > 0)
0651 herun_->Fill(rnnum_, (double)(count2HE) / countHE);
0652 }
0653 #ifdef EDM_ML_DEBUG
0654 edm::LogVerbatim("RecAnalyzerMinbias") << "HBHE " << count2 << ":" << count << ":" << (double)(count2) / count
0655 << "\t HB " << count2HB << ":" << countHB << ":"
0656 << (double)(count2HB) / countHB << "\t HE " << count2HE << ":" << countHE
0657 << ":" << (double)(count2HE) / countHE;
0658 #endif
0659 int countHF(0), count2HF(0);
0660
0661 for (HFRecHitCollection::const_iterator hfItr = HithfMB.begin(); hfItr != HithfMB.end(); hfItr++) {
0662
0663 DetId mydetid = hfItr->id().rawId();
0664 double icalconst(1.);
0665 if (theRecalib_) {
0666 std::map<DetId, double>::iterator itr = corrFactor_.find(mydetid);
0667 if (itr != corrFactor_.end())
0668 icalconst = itr->second;
0669 }
0670 HFRecHit aHit(hfItr->id(), hfItr->energy() * icalconst, hfItr->time());
0671
0672 double energyhit = aHit.energy();
0673 DetId id = (*hfItr).detid();
0674 HcalDetId hid = HcalDetId(id);
0675 ++countHF;
0676 if (fill) {
0677 for (unsigned int i = 0; i < hcalID_.size(); i++) {
0678 if (hcalID_[i] == id.rawId()) {
0679 histo_[i]->Fill(energyhit);
0680 break;
0681 }
0682 }
0683 if (fillHist_) {
0684 std::map<HcalDetId, TH1D*>::iterator itr1 = histHC_.find(hid);
0685 if (itr1 != histHC_.end())
0686 itr1->second->Fill(energyhit);
0687 }
0688 h_[hid.subdet() - 1]->Fill(energyhit);
0689 if (energyhit > eMin_) {
0690 hf_->Fill(hid.ieta(), hid.iphi());
0691 ++count2HF;
0692 }
0693 }
0694
0695
0696
0697
0698 if (!fillHist_) {
0699 if (((Noise_ || runNZS_) && fabs(energyhit) <= 40.) || (energyhit >= eLowHF_ && energyhit <= eHighHF_)) {
0700 std::map<std::pair<int, HcalDetId>, myInfo>::iterator itr1 =
0701 myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0702 if (itr1 == myMap_.end()) {
0703 myInfo info;
0704 myMap_[std::pair<int, HcalDetId>(algoBit, hid)] = info;
0705 itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0706 }
0707 itr1->second.theMB0 += weight;
0708 itr1->second.theMB1 += (weight * energyhit);
0709 itr1->second.theMB2 += (weight * energyhit * energyhit);
0710 itr1->second.theMB3 += (weight * energyhit * energyhit * energyhit);
0711 itr1->second.theMB4 += (weight * energyhit * energyhit * energyhit * energyhit);
0712 itr1->second.runcheck = rnnum_;
0713 }
0714 }
0715 }
0716 if (fill && countHF > 0)
0717 hfrun_->Fill(rnnum_, (double)(count2HF) / countHF);
0718 #ifdef EDM_ML_DEBUG
0719 if (count)
0720 edm::LogVerbatim("RecAnalyzerMinbias")
0721 << "HF " << count2HF << ":" << countHF << ":" << (double)(count2HF) / countHF;
0722 #endif
0723 }
0724
0725
0726 #include "FWCore/Framework/interface/MakerMacros.h"
0727 DEFINE_FWK_MODULE(RecAnalyzerMinbias);