HcalCalibAlgos/plugins/SimAnalyzerMinbias.cc

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 "FWCore/Framework/interface/Frameworkfwd.h"
0012 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0013 #include "FWCore/Framework/interface/Event.h"
0014 #include "FWCore/Framework/interface/EventSetup.h"
0015 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0016 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0017 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0018 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0019 #include "FWCore/ServiceRegistry/interface/Service.h"
0020 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0021
0022 #include "CondFormats/HcalObjects/interface/HcalRespCorrs.h"
0023 #include "CondFormats/DataRecord/interface/HcalRespCorrsRcd.h"
0024 #include "CondFormats/DataRecord/interface/L1GtTriggerMenuRcd.h"
0025 #include "CondFormats/L1TObjects/interface/L1GtTriggerMenu.h"
0026
0027 #include "DataFormats/DetId/interface/DetId.h"
0028 #include "DataFormats/HcalDetId/interface/HcalSubdetector.h"
0029 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0030 #include "SimDataFormats/CaloHit/interface/PCaloHitContainer.h"
0031 #include "SimDataFormats/GeneratorProducts/interface/HepMCProduct.h"
0032
0033 #include "EventFilter/HcalRawToDigi/interface/HcalDCCHeader.h"
0034
0035 #include "TFile.h"
0036 #include "TTree.h"
0037
0038 // class declaration
0039 class SimAnalyzerMinbias : public edm::one::EDAnalyzer<edm::one::WatchRuns, edm::one::SharedResources> {
0040 public:
0041   explicit SimAnalyzerMinbias(const edm::ParameterSet&);
0042   static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0043   void analyze(const edm::Event&, const edm::EventSetup&) override;
0044   void beginJob() override;
0045   void endJob() override;
0046   void beginRun(const edm::Run& r, const edm::EventSetup& iSetup) override{};
0047   void endRun(const edm::Run& r, const edm::EventSetup& iSetup) override{};
0048
0049 private:
0050   // ----------member data ---------------------------
0051   const double timeCut_;
0052   TTree* myTree_;
0053
0054   // Root tree members
0055   int mydet, mysubd, depth, iphi, ieta, cells;
0056   float mom0_MB, mom1_MB, mom2_MB, mom3_MB, mom4_MB;
0057   struct myInfo {
0058     double theMB0, theMB1, theMB2, theMB3, theMB4;
0059     myInfo() { theMB0 = theMB1 = theMB2 = theMB3 = theMB4 = 0; }
0060   };
0061   std::map<HcalDetId, myInfo> myMap_;
0062   const edm::EDGetTokenT<edm::HepMCProduct> tok_evt_;
0063   const edm::EDGetTokenT<edm::PCaloHitContainer> tok_hcal_;
0064 };
0065
0066 // constructors and destructor
0067
0068 SimAnalyzerMinbias::SimAnalyzerMinbias(const edm::ParameterSet& iConfig)
0069     : timeCut_(iConfig.getUntrackedParameter<double>("TimeCut", 500)),
0070       tok_evt_(consumes<edm::HepMCProduct>(edm::InputTag("generator"))),
0071       tok_hcal_(consumes<edm::PCaloHitContainer>(edm::InputTag("g4SimHits", "HcalHits"))) {
0072   usesResource(TFileService::kSharedResource);
0073   edm::LogVerbatim("AnalyzerMB") << "Use Time cut of " << timeCut_ << " ns";
0074 }
0075
0076 void SimAnalyzerMinbias::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0077   edm::ParameterSetDescription desc;
0078   desc.addUntracked<double>("TimeCut", 500);
0079   descriptions.add("simAnalyzerMinbias", desc);
0080 }
0081
0082 void SimAnalyzerMinbias::beginJob() {
0083   edm::Service<TFileService> fs;
0084   myTree_ = fs->make<TTree>("SimJet", "SimJet Tree");
0085   myTree_->Branch("mydet", &mydet, "mydet/I");
0086   myTree_->Branch("mysubd", &mysubd, "mysubd/I");
0087   myTree_->Branch("cells", &cells, "cells");
0088   myTree_->Branch("depth", &depth, "depth/I");
0089   myTree_->Branch("ieta", &ieta, "ieta/I");
0090   myTree_->Branch("iphi", &iphi, "iphi/I");
0091   myTree_->Branch("mom0_MB", &mom0_MB, "mom0_MB/F");
0092   myTree_->Branch("mom1_MB", &mom1_MB, "mom1_MB/F");
0093   myTree_->Branch("mom2_MB", &mom2_MB, "mom2_MB/F");
0094   myTree_->Branch("mom3_MB", &mom3_MB, "mom3_MB/F");
0095   myTree_->Branch("mom4_MB", &mom4_MB, "mom4_MB/F");
0096
0097   myMap_.clear();
0098 }
0099
0100 //  EndJob
0101 //
0102 void SimAnalyzerMinbias::endJob() {
0103   cells = 0;
0104   for (std::map<HcalDetId, myInfo>::const_iterator itr = myMap_.begin(); itr != myMap_.end(); ++itr) {
0105     mysubd = itr->first.subdet();
0106     depth = itr->first.depth();
0107     iphi = itr->first.iphi();
0108     ieta = itr->first.ieta();
0109     myInfo info = itr->second;
0110     if (info.theMB0 > 0) {
0111       mom0_MB = info.theMB0;
0112       mom1_MB = info.theMB1;
0113       mom2_MB = info.theMB2;
0114       mom3_MB = info.theMB3;
0115       mom4_MB = info.theMB4;
0116       cells++;
0117
0118       edm::LogVerbatim("SimAnalyzerMinbias")
0119           << " Result=  " << mysubd << " " << ieta << " " << iphi << " mom0  " << mom0_MB << " mom1 " << mom1_MB
0120           << " mom2 " << mom2_MB << " mom3 " << mom3_MB << " mom4 " << mom4_MB;
0121       myTree_->Fill();
0122     }
0123   }
0124   edm::LogVerbatim("SimAnalyzerMinbias") << "cells " << cells;
0125 }
0126
0127 //
0128 // member functions
0129 //
0130
0131 // ------------ method called to produce the data  ------------
0132
0133 void SimAnalyzerMinbias::analyze(const edm::Event& iEvent, const edm::EventSetup&) {
0134   edm::LogVerbatim("SimAnalyzerMinbias") << " Start SimAnalyzerMinbias::analyze " << iEvent.id().run() << ":"
0135                                          << iEvent.id().event();
0136
0137   const edm::Handle<edm::HepMCProduct>& evtMC = iEvent.getHandle(tok_evt_);
0138   if (!evtMC.isValid()) {
0139     edm::LogWarning("SimAnalyzerMinbias") << "no HepMCProduct found";
0140   } else {
0141     const HepMC::GenEvent* myGenEvent = evtMC->GetEvent();
0142     edm::LogVerbatim("SimAnalyzerMinbias") << "Event with " << myGenEvent->particles_size() << " particles + "
0143                                            << myGenEvent->vertices_size() << " vertices";
0144   }
0145
0146   const edm::Handle<edm::PCaloHitContainer>& hcalHits = iEvent.getHandle(tok_hcal_);
0147   if (!hcalHits.isValid()) {
0148     edm::LogWarning("SimAnalyzerMinbias") << "Error! can't get HcalHits product!";
0149     return;
0150   }
0151
0152   const edm::PCaloHitContainer* HitHcal = hcalHits.product();
0153   std::map<HcalDetId, double> hitMap;
0154   for (std::vector<PCaloHit>::const_iterator hcalItr = HitHcal->begin(); hcalItr != HitHcal->end(); ++hcalItr) {
0155     double time = hcalItr->time();
0156     if (time < timeCut_) {
0157       double energyhit = hcalItr->energy();
0158       HcalDetId hid = HcalDetId(hcalItr->id());
0159       std::map<HcalDetId, double>::iterator itr1 = hitMap.find(hid);
0160       if (itr1 == hitMap.end()) {
0161         hitMap[hid] = 0;
0162         itr1 = hitMap.find(hid);
0163       }
0164       itr1->second += energyhit;
0165     }
0166   }
0167   edm::LogVerbatim("SimAnalyzerMinbias") << "extract information of " << hitMap.size() << " towers from "
0168                                          << HitHcal->size() << " hits";
0169
0170   for (std::map<HcalDetId, double>::const_iterator hcalItr = hitMap.begin(); hcalItr != hitMap.end(); ++hcalItr) {
0171     HcalDetId hid = hcalItr->first;
0172     double energyhit = hcalItr->second;
0173     std::map<HcalDetId, myInfo>::iterator itr1 = myMap_.find(hid);
0174     if (itr1 == myMap_.end()) {
0175       myInfo info;
0176       myMap_[hid] = info;
0177       itr1 = myMap_.find(hid);
0178     }
0179     itr1->second.theMB0++;
0180     itr1->second.theMB1 += energyhit;
0181     itr1->second.theMB2 += (energyhit * energyhit);
0182     itr1->second.theMB3 += (energyhit * energyhit * energyhit);
0183     itr1->second.theMB4 += (energyhit * energyhit * energyhit * energyhit);
0184     edm::LogVerbatim("SimAnalyzerMinbias") << "ID " << hid << " with energy " << energyhit;
0185   }
0186 }
0187
0188 //define this as a plug-in
0189 #include "FWCore/Framework/interface/MakerMacros.h"
0190 DEFINE_FWK_MODULE(SimAnalyzerMinbias);