File indexing completed on 2024-09-07 04:34:59
0001
0002 #include <memory>
0003 #include <string>
0004 #include <iostream>
0005 #include <fstream>
0006 #include <sstream>
0007 #include <vector>
0008 #include <map>
0009
0010
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
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
0051 const double timeCut_;
0052 TTree* myTree_;
0053
0054
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
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
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
0129
0130
0131
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
0189 #include "FWCore/Framework/interface/MakerMacros.h"
0190 DEFINE_FWK_MODULE(SimAnalyzerMinbias);