Back to home page

Project CMSSW displayed by LXR
 
 

    

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

0001 // system include files
0002 #include <memory>
0003 #include <string>
0004 #include <iostream>
0005 #include <fstream>
0006 #include <sstream>
0007 #include <vector>
0008 #include <map>
0009 
0010 // user include files
0011 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0012 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0013 #include "FWCore/Common/interface/TriggerNames.h"
0014 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0015 #include "FWCore/Framework/interface/Event.h"
0016 #include "FWCore/Framework/interface/EventSetup.h"
0017 #include "FWCore/Framework/interface/Frameworkfwd.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0021 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0022 #include "FWCore/ServiceRegistry/interface/Service.h"
0023 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0024 #include "DataFormats/Common/interface/TriggerResults.h"
0025 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0026 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0027 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0028 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
0029 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
0030 #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
0031 
0032 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0033 
0034 #include "TH1D.h"
0035 #include "TTree.h"
0036 
0037 //#define EDM_ML_DEBUG
0038 
0039 // class declaration
0040 class RecAnalyzerHF : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0041 public:
0042   explicit RecAnalyzerHF(const edm::ParameterSet&);
0043   ~RecAnalyzerHF() override = default;
0044 
0045   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0046 
0047 private:
0048   void analyze(edm::Event const&, edm::EventSetup const&) override;
0049   void beginJob() override;
0050   void endJob() override;
0051 
0052 private:
0053   void analyzeHcal(const HFPreRecHitCollection&, int, bool);
0054 
0055   // ----------member data ---------------------------
0056   bool nzs_, noise_, ratio_, ignoreL1_, fillTree_;
0057   double eLowHF_, eHighHF_;
0058   std::vector<int> trigbit_;
0059   std::vector<unsigned int> hcalID_;
0060   TTree* myTree_;
0061   TH1D* hist_[2];
0062   std::vector<std::pair<TH1D*, TH1D*>> histo_;
0063   double rnnum_;
0064   struct myInfo {
0065     double kount, f11, f12, f13, f14, f21, f22, f23, f24, runcheck;
0066     myInfo() { kount = f11 = f12 = f13 = f14 = f21 = f22 = f23 = f24 = runcheck = 0; }
0067   };
0068   // Root tree members
0069   double rnnumber;
0070   int mysubd, depth, iphi, ieta, cells, trigbit;
0071   float mom0_F1, mom1_F1, mom2_F1, mom3_F1, mom4_F1;
0072   float mom0_F2, mom1_F2, mom2_F2, mom3_F2, mom4_F2;
0073   std::map<std::pair<int, HcalDetId>, myInfo> myMap_;
0074   edm::EDGetTokenT<HFPreRecHitCollection> tok_hfreco_;
0075   edm::EDGetTokenT<L1GlobalTriggerObjectMapRecord> tok_hltL1GtMap_;
0076 };
0077 
0078 // constructors and destructor
0079 RecAnalyzerHF::RecAnalyzerHF(const edm::ParameterSet& iConfig)
0080     : nzs_(iConfig.getParameter<bool>("RunNZS")),
0081       noise_(iConfig.getParameter<bool>("Noise")),
0082       ratio_(iConfig.getParameter<bool>("Ratio")),
0083       ignoreL1_(iConfig.getUntrackedParameter<bool>("IgnoreL1", false)),
0084       fillTree_(iConfig.getUntrackedParameter<bool>("FillTree", true)),
0085       eLowHF_(iConfig.getParameter<double>("ELowHF")),
0086       eHighHF_(iConfig.getParameter<double>("EHighHF")),
0087       trigbit_(iConfig.getUntrackedParameter<std::vector<int>>("TriggerBits")),
0088       tok_hfreco_(consumes<HFPreRecHitCollection>(iConfig.getParameter<edm::InputTag>("hfInput"))),
0089       tok_hltL1GtMap_(consumes<L1GlobalTriggerObjectMapRecord>(edm::InputTag("hltL1GtObjectMap"))) {
0090   usesResource("TFileService");
0091   // get constants for the ID
0092   std::vector<int> ieta = iConfig.getUntrackedParameter<std::vector<int>>("HcalIeta");
0093   std::vector<int> iphi = iConfig.getUntrackedParameter<std::vector<int>>("HcalIphi");
0094   std::vector<int> depth = iConfig.getUntrackedParameter<std::vector<int>>("HcalDepth");
0095 
0096   edm::LogVerbatim("RecAnalyzerHF") << " Flags (IgnoreL1): " << ignoreL1_ << " (NZS) " << nzs_ << " (Noise) " << noise_
0097                                     << " (Ratio) " << ratio_;
0098   edm::LogVerbatim("RecAnalyzerHF") << "Thresholds for HF " << eLowHF_ << ":" << eHighHF_;
0099   for (unsigned int k = 0; k < ieta.size(); ++k) {
0100     if (std::abs(ieta[k]) >= 29) {
0101       unsigned int id = (HcalDetId(HcalForward, ieta[k], iphi[k], depth[k])).rawId();
0102       hcalID_.push_back(id);
0103       edm::LogVerbatim("RecAnalyzerHF") << "DetId[" << k << "] " << HcalDetId(id);
0104     }
0105   }
0106   edm::LogVerbatim("RecAnalyzerHF") << "Select on " << trigbit_.size() << " L1 Trigger selection";
0107   unsigned int k(0);
0108   for (auto trig : trigbit_) {
0109     edm::LogVerbatim("RecAnalyzerHF") << "Bit[" << k << "] " << trig;
0110     ++k;
0111   }
0112 }
0113 
0114 // ------------ method fills 'descriptions' with the allowed parameters for the module  ------------
0115 void RecAnalyzerHF::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0116   edm::ParameterSetDescription desc;
0117   desc.add<bool>("RunNZS", true);
0118   desc.add<bool>("Noise", false);
0119   desc.add<bool>("Ratio", false);
0120   desc.add<double>("ELowHF", 10);
0121   desc.add<double>("EHighHF", 150);
0122   std::vector<int> idummy;
0123   desc.addUntracked<std::vector<int>>("TriggerBits", idummy);
0124   desc.addUntracked<bool>("IgnoreL1", false);
0125   desc.addUntracked<bool>("FillHisto", false);
0126   desc.addUntracked<std::vector<int>>("HcalIeta", idummy);
0127   desc.addUntracked<std::vector<int>>("HcalIphi", idummy);
0128   desc.addUntracked<std::vector<int>>("HcalDepth", idummy);
0129   desc.add<edm::InputTag>("hfInput", edm::InputTag("hfprereco"));
0130   descriptions.add("recAnalyzerHF", desc);
0131 }
0132 
0133 void RecAnalyzerHF::beginJob() {
0134   edm::Service<TFileService> fs;
0135   char name[700], title[700];
0136   double xmin(-10.0), xmax(90.0);
0137   if (ratio_) {
0138     xmin = -5.0;
0139     xmax = 5.0;
0140   }
0141   for (int i = 0; i < 2; ++i) {
0142     sprintf(name, "HF%d", i);
0143     sprintf(title, "The metric F%d for HF", i + 1);
0144     hist_[i] = fs->make<TH1D>(name, title, 50, xmin, xmax);
0145   }
0146 
0147   for (const auto& id : hcalID_) {
0148     HcalDetId hid = HcalDetId(id);
0149     TH1D *h1(nullptr), *h2(nullptr);
0150     for (int i = 0; i < 2; ++i) {
0151       sprintf(name, "HF%d%d_%d_%d", i, hid.ieta(), hid.iphi(), hid.depth());
0152       sprintf(title, "The metric F%d for HF i#eta %d i#phi %d depth %d", i + 1, hid.ieta(), hid.iphi(), hid.depth());
0153       if (i == 0)
0154         h1 = fs->make<TH1D>(name, title, 50, xmin, xmax);
0155       else
0156         h2 = fs->make<TH1D>(name, title, 50, xmin, xmax);
0157     }
0158     histo_.push_back(std::pair<TH1D*, TH1D*>(h1, h2));
0159   };
0160 
0161   if (fillTree_) {
0162     myTree_ = fs->make<TTree>("RecJet", "RecJet Tree");
0163     myTree_->Branch("cells", &cells, "cells/I");
0164     myTree_->Branch("mysubd", &mysubd, "mysubd/I");
0165     myTree_->Branch("depth", &depth, "depth/I");
0166     myTree_->Branch("ieta", &ieta, "ieta/I");
0167     myTree_->Branch("iphi", &iphi, "iphi/I");
0168     myTree_->Branch("mom0_F1", &mom0_F1, "mom0_F1/F");
0169     myTree_->Branch("mom1_F1", &mom1_F1, "mom1_F1/F");
0170     myTree_->Branch("mom2_F1", &mom2_F1, "mom2_F1/F");
0171     myTree_->Branch("mom3_F1", &mom3_F1, "mom3_F1/F");
0172     myTree_->Branch("mom4_F1", &mom4_F1, "mom4_F1/F");
0173     myTree_->Branch("mom0_F2", &mom0_F2, "mom0_F2/F");
0174     myTree_->Branch("mom1_F2", &mom1_F2, "mom1_F2/F");
0175     myTree_->Branch("mom2_F2", &mom2_F2, "mom2_F2/F");
0176     myTree_->Branch("mom3_F2", &mom3_F2, "mom3_F2/F");
0177     myTree_->Branch("mom4_F2", &mom4_F2, "mom4_F2/F");
0178     myTree_->Branch("trigbit", &trigbit, "trigbit/I");
0179     myTree_->Branch("rnnumber", &rnnumber, "rnnumber/D");
0180   }
0181 
0182   myMap_.clear();
0183 }
0184 
0185 void RecAnalyzerHF::endJob() {
0186   if (fillTree_) {
0187     cells = 0;
0188     for (const auto& itr : myMap_) {
0189       edm::LogVerbatim("RecAnalyzerHF") << "Fired trigger bit number " << itr.first.first;
0190       myInfo info = itr.second;
0191       if (info.kount > 0) {
0192         mom0_F1 = info.kount;
0193         mom1_F1 = info.f11;
0194         mom2_F1 = info.f12;
0195         mom3_F1 = info.f13;
0196         mom4_F1 = info.f14;
0197         mom0_F2 = info.kount;
0198         mom1_F2 = info.f21;
0199         mom2_F2 = info.f22;
0200         mom3_F2 = info.f23;
0201         mom4_F2 = info.f24;
0202         rnnumber = info.runcheck;
0203         trigbit = itr.first.first;
0204         mysubd = itr.first.second.subdet();
0205         depth = itr.first.second.depth();
0206         iphi = itr.first.second.iphi();
0207         ieta = itr.first.second.ieta();
0208         edm::LogVerbatim("RecAnalyzerHF")
0209             << " Result=  " << trigbit << " " << mysubd << " " << ieta << " " << iphi << " F1:mom0  " << mom0_F1
0210             << " mom1 " << mom1_F1 << " mom2 " << mom2_F1 << " mom3 " << mom3_F1 << " mom4 " << mom4_F1 << " F2:mom0 "
0211             << mom0_F2 << " mom1 " << mom1_F2 << " mom2 " << mom2_F2 << " mom3 " << mom3_F2 << " mom4 " << mom4_F2;
0212         myTree_->Fill();
0213         cells++;
0214       }
0215     }
0216     edm::LogVerbatim("RecAnalyzerHF") << "cells"
0217                                       << " " << cells;
0218   }
0219 #ifdef EDM_ML_DEBUG
0220   edm::LogVerbatim("RecAnalyzerHF") << "Exiting from RecAnalyzerHF::endjob";
0221 #endif
0222 }
0223 
0224 //
0225 // member functions
0226 //
0227 // ------------ method called to produce the data  ------------
0228 
0229 void RecAnalyzerHF::analyze(const edm::Event& iEvent, const edm::EventSetup&) {
0230   rnnum_ = (double)iEvent.run();
0231 
0232   const edm::Handle<HFPreRecHitCollection>& hf = iEvent.getHandle(tok_hfreco_);
0233   if (!hf.isValid()) {
0234     edm::LogWarning("RecAnalyzerHF") << "HcalCalibAlgos: Error! can't get hf product!";
0235     return;
0236   }
0237   const HFPreRecHitCollection Hithf = *(hf.product());
0238   edm::LogVerbatim("RecAnalyzerHF") << "HF MB size of collection " << Hithf.size();
0239   if (Hithf.size() < 1700 && nzs_) {
0240     edm::LogWarning("RecAnalyzerHF") << "HF problem " << rnnum_ << " size " << Hithf.size();
0241   }
0242 
0243   bool select(false);
0244   if (!trigbit_.empty()) {
0245     const edm::Handle<L1GlobalTriggerObjectMapRecord>& gtObjectMapRecord = iEvent.getHandle(tok_hltL1GtMap_);
0246     if (gtObjectMapRecord.isValid()) {
0247       const std::vector<L1GlobalTriggerObjectMap>& objMapVec = gtObjectMapRecord->gtObjectMap();
0248       for (std::vector<L1GlobalTriggerObjectMap>::const_iterator itMap = objMapVec.begin(); itMap != objMapVec.end();
0249            ++itMap) {
0250         bool resultGt = (*itMap).algoGtlResult();
0251         if (resultGt) {
0252           int algoBit = (*itMap).algoBitNumber();
0253           if (std::find(trigbit_.begin(), trigbit_.end(), algoBit) != trigbit_.end()) {
0254             select = true;
0255             break;
0256           }
0257         }
0258       }
0259     }
0260   }
0261 
0262   //event weight for FLAT sample and PU information
0263 
0264   if (ignoreL1_ || (!trigbit_.empty() && select)) {
0265     analyzeHcal(Hithf, 1, true);
0266   } else if ((!ignoreL1_) && (trigbit_.empty())) {
0267     const edm::Handle<L1GlobalTriggerObjectMapRecord>& gtObjectMapRecord = iEvent.getHandle(tok_hltL1GtMap_);
0268     if (gtObjectMapRecord.isValid()) {
0269       const std::vector<L1GlobalTriggerObjectMap>& objMapVec = gtObjectMapRecord->gtObjectMap();
0270       bool ok(false);
0271       for (std::vector<L1GlobalTriggerObjectMap>::const_iterator itMap = objMapVec.begin(); itMap != objMapVec.end();
0272            ++itMap) {
0273         bool resultGt = (*itMap).algoGtlResult();
0274         if (resultGt) {
0275           int algoBit = (*itMap).algoBitNumber();
0276           analyzeHcal(Hithf, algoBit, (!ok));
0277           ok = true;
0278         }
0279       }
0280       if (!ok) {
0281         edm::LogVerbatim("RecAnalyzerHF") << "No passed L1 Trigger found";
0282       }
0283     }
0284   }
0285 }
0286 
0287 void RecAnalyzerHF::analyzeHcal(const HFPreRecHitCollection& Hithf, int algoBit, bool fill) {
0288 #ifdef EDM_ML_DEBUG
0289   edm::LogVerbatim("RecAnalyzerHF") << "Enter analyzeHcal for bit " << algoBit << " Fill " << fill
0290                                     << " Collection size " << Hithf.size();
0291 #endif
0292   // Signal part for HF
0293   for (const auto& hfItr : Hithf) {
0294     HcalDetId hid = hfItr.id();
0295     double e0 = (hfItr.getHFQIE10Info(0) == nullptr) ? 0 : hfItr.getHFQIE10Info(0)->energy();
0296     double e1 = (hfItr.getHFQIE10Info(1) == nullptr) ? 0 : hfItr.getHFQIE10Info(1)->energy();
0297     double energy = e0 + e1;
0298     if (std::abs(energy) < 1e-6)
0299       energy = (energy > 0) ? 1e-6 : -1e-6;
0300     double f1(e0), f2(e1);
0301     if (ratio_) {
0302       f1 /= energy;
0303       f2 /= energy;
0304     }
0305 #ifdef EDM_ML_DEBUG
0306     edm::LogVerbatim("RecAnalyzerHF") << hid << " E " << e0 << ":" << e1 << " F " << f1 << ":" << f2;
0307 #endif
0308     if (fill) {
0309       for (unsigned int i = 0; i < hcalID_.size(); i++) {
0310         if (hcalID_[i] == hid.rawId()) {
0311           histo_[i].first->Fill(f1);
0312           histo_[i].second->Fill(f2);
0313           break;
0314         }
0315       }
0316       hist_[0]->Fill(f1);
0317       hist_[1]->Fill(f2);
0318     }
0319 
0320     //
0321     // Remove PMT hits
0322     //
0323     if (((noise_ || nzs_) && fabs(energy) <= 40.) || (energy >= eLowHF_ && energy <= eHighHF_)) {
0324       std::map<std::pair<int, HcalDetId>, myInfo>::iterator itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0325       if (itr1 == myMap_.end()) {
0326         myInfo info;
0327         myMap_[std::pair<int, HcalDetId>(algoBit, hid)] = info;
0328         itr1 = myMap_.find(std::pair<int, HcalDetId>(algoBit, hid));
0329       }
0330       itr1->second.kount += 1.0;
0331       itr1->second.f11 += (f1);
0332       itr1->second.f12 += (f1 * f1);
0333       itr1->second.f13 += (f1 * f1 * f1);
0334       itr1->second.f14 += (f1 * f1 * f1 * f1);
0335       itr1->second.f21 += (f2);
0336       itr1->second.f22 += (f2 * f2);
0337       itr1->second.f23 += (f2 * f2 * f2);
0338       itr1->second.f24 += (f2 * f2 * f2 * f2);
0339       itr1->second.runcheck = rnnum_;
0340     }
0341   }
0342 }
0343 
0344 //define this as a plug-in
0345 #include "FWCore/Framework/interface/MakerMacros.h"
0346 DEFINE_FWK_MODULE(RecAnalyzerHF);
[ 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!