Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 11:59:03 

0001 // -*- C++ -*-
0002 //
0003 // Package:    HcalCalibAlgos
0004 // Class:      RecAnalyzerMinbias
0005 //
0006 /**\class RecAnalyzerMinbias RecAnalyzerMinbias.cc Calibration/HcalCalibAlgos/test/RecAnalyzerMinbias.cc
0007 
0008  Description: Performs phi-symmetry studies of HB/HE/HF channels
0009 
0010  Implementation:
0011      <Notes on implementation>
0012 */
0013 //
0014 // Original Author:  Sunanda Banerjee
0015 //         Created:  Thu Mar  4 18:52:02 CST 2012
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 #include <string>
0022 #include <iostream>
0023 #include <fstream>
0024 #include <sstream>
0025 #include <vector>
0026 #include <map>
0027 
0028 // user include files
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 //#define EDM_ML_DEBUG
0065 
0066 // class declaration
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   // ----------member data ---------------------------
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   // Root tree members
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 // constructors and destructor
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   // get constants for DetID's
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   // Read correction factors
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   // Suitable cutoff to remove fluctuation of pedestal
0223   desc.addUntracked<double>("eMin", 2.0);
0224   // The following run range is suited to study 2017 commissioning period
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 //  EndJob
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       // For HB
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       // For HE
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       // For HF
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 // member functions
0415 //
0416 // ------------ method called to produce the data  ------------
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   //event weight for FLAT sample and PU information
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   // Signal part for HB HE
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     // Recalibration of energy
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   }  // HBHE_MB
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   // Signal part for HF
0661   for (HFRecHitCollection::const_iterator hfItr = HithfMB.begin(); hfItr != HithfMB.end(); hfItr++) {
0662     // Recalibration of energy
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     // Remove PMT hits
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 //define this as a plug-in
0726 #include "FWCore/Framework/interface/MakerMacros.h"
0727 DEFINE_FWK_MODULE(RecAnalyzerMinbias);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!