File indexing completed on 2023-03-17 11:10:45
0001
0002 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0004 #include "FWCore/ServiceRegistry/interface/Service.h"
0005 #include "JetMETCorrections/MinBias/interface/MinBias.h"
0006 #include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"
0007
0008 namespace cms {
0009 MinBias::MinBias(const edm::ParameterSet& iConfig)
0010 : caloGeometryESToken_(esConsumes<edm::Transition::BeginRun>()), geo_(nullptr) {
0011
0012 hbheLabel_ = iConfig.getParameter<std::string>("hbheInput");
0013 hoLabel_ = iConfig.getParameter<std::string>("hoInput");
0014 hfLabel_ = iConfig.getParameter<std::string>("hfInput");
0015 hbheToken_ = mayConsume<HBHERecHitCollection>(edm::InputTag(hbheLabel_));
0016 hoToken_ = mayConsume<HORecHitCollection>(edm::InputTag(hoLabel_));
0017 hfToken_ = mayConsume<HFRecHitCollection>(edm::InputTag(hfLabel_));
0018 allowMissingInputs_ = iConfig.getUntrackedParameter<bool>("AllowMissingInputs", false);
0019 }
0020
0021 void MinBias::beginJob() {
0022 edm::Service<TFileService> fs;
0023 myTree_ = fs->make<TTree>("RecJet", "RecJet Tree");
0024 myTree_->Branch("mydet", &mydet, "mydet/I");
0025 myTree_->Branch("mysubd", &mysubd, "mysubd/I");
0026 myTree_->Branch("depth", &depth, "depth/I");
0027 myTree_->Branch("ieta", &ieta, "ieta/I");
0028 myTree_->Branch("iphi", &iphi, "iphi/I");
0029 myTree_->Branch("eta", &eta, "eta/F");
0030 myTree_->Branch("phi", &phi, "phi/F");
0031 myTree_->Branch("mom1", &mom1, "mom1/F");
0032 myTree_->Branch("mom2", &mom2, "mom2/F");
0033 myTree_->Branch("mom3", &mom3, "mom3/F");
0034 myTree_->Branch("mom4", &mom4, "mom4/F");
0035 }
0036
0037 void MinBias::beginRun(edm::Run const&, edm::EventSetup const& iSetup) {
0038 geo_ = &iSetup.getData(caloGeometryESToken_);
0039 std::vector<DetId> did = geo_->getValidDetIds();
0040
0041 for (auto const& id : did) {
0042 if ((id).det() == DetId::Hcal) {
0043 theFillDetMap0_[id] = 0.;
0044 theFillDetMap1_[id] = 0.;
0045 theFillDetMap2_[id] = 0.;
0046 theFillDetMap3_[id] = 0.;
0047 theFillDetMap4_[id] = 0.;
0048 }
0049 }
0050 }
0051
0052 void MinBias::endRun(edm::Run const&, edm::EventSetup const& iSetup) {}
0053
0054 void MinBias::endJob() {
0055 const HcalGeometry* hgeo = static_cast<const HcalGeometry*>(geo_->getSubdetectorGeometry(DetId::Hcal, 1));
0056 const std::vector<DetId>& did = hgeo->getValidDetIds();
0057 int i = 0;
0058 for (const auto& id : did) {
0059
0060 GlobalPoint pos = hgeo->getPosition(id);
0061 mydet = (int)(id.det());
0062 mysubd = (id.subdetId());
0063 depth = HcalDetId(id).depth();
0064 ieta = HcalDetId(id).ieta();
0065 iphi = HcalDetId(id).iphi();
0066 phi = pos.phi();
0067 eta = pos.eta();
0068 if (theFillDetMap0_[id] > 0.) {
0069 mom1 = theFillDetMap1_[id] / theFillDetMap0_[id];
0070 mom2 = theFillDetMap2_[id] / theFillDetMap0_[id] - (mom1 * mom1);
0071 mom3 = theFillDetMap3_[id] / theFillDetMap0_[id] - 3. * mom1 * theFillDetMap2_[id] / theFillDetMap0_[id] +
0072 2. * pow(mom2, 3);
0073 mom4 = (theFillDetMap4_[id] - 4. * mom1 * theFillDetMap3_[id] + 6. * pow(mom1, 2) * theFillDetMap2_[id]) /
0074 theFillDetMap0_[id] -
0075 3. * pow(mom1, 4);
0076
0077 } else {
0078 mom1 = 0.;
0079 mom2 = 0.;
0080 mom3 = 0.;
0081 mom4 = 0.;
0082 }
0083 edm::LogWarning("MinBias") << " Detector " << id.rawId() << " mydet " << mydet << " " << mysubd << " " << depth
0084 << " " << ieta << " " << iphi << " " << pos.eta() << " " << pos.phi();
0085 edm::LogWarning("MinBias") << " Energy " << mom1 << " " << mom2 << std::endl;
0086 myTree_->Fill();
0087 i++;
0088
0089 }
0090 edm::LogWarning("MinBias") << " The number of CaloDet records " << did.size();
0091 edm::LogWarning("MinBias") << " The number of Hcal records " << i << std::endl;
0092
0093
0094
0095
0096
0097
0098
0099
0100
0101 }
0102
0103
0104
0105
0106
0107
0108 void MinBias::analyze(const edm::Event& iEvent, const edm::EventSetup&) {
0109 if (!hbheLabel_.empty()) {
0110 edm::Handle<HBHERecHitCollection> hbhe;
0111 iEvent.getByToken(hbheToken_, hbhe);
0112 if (!hbhe.isValid()) {
0113
0114 if (!allowMissingInputs_) {
0115 *hbhe;
0116 }
0117 } else {
0118 for (auto const& hbheItr : (HBHERecHitCollection)(*hbhe)) {
0119 DetId id = (hbheItr).detid();
0120 if (hbheItr.energy() > 0.)
0121 edm::LogWarning("MinBias") << " Energy = " << hbheItr.energy();
0122 theFillDetMap0_[id] += 1.;
0123 theFillDetMap1_[id] += hbheItr.energy();
0124 theFillDetMap2_[id] += pow(hbheItr.energy(), 2);
0125 theFillDetMap3_[id] += pow(hbheItr.energy(), 3);
0126 theFillDetMap4_[id] += pow(hbheItr.energy(), 4);
0127 }
0128 }
0129 }
0130
0131 if (!hoLabel_.empty()) {
0132 edm::Handle<HORecHitCollection> ho;
0133 iEvent.getByToken(hoToken_, ho);
0134 if (!ho.isValid()) {
0135
0136 if (!allowMissingInputs_) {
0137 *ho;
0138 }
0139 } else {
0140 for (auto const& hoItr : (HORecHitCollection)(*ho)) {
0141 DetId id = hoItr.detid();
0142 theFillDetMap0_[id] += 1.;
0143 theFillDetMap1_[id] += hoItr.energy();
0144 theFillDetMap2_[id] += pow(hoItr.energy(), 2);
0145 theFillDetMap3_[id] += pow(hoItr.energy(), 3);
0146 theFillDetMap4_[id] += pow(hoItr.energy(), 4);
0147 }
0148 }
0149 }
0150
0151 if (!hfLabel_.empty()) {
0152 edm::Handle<HFRecHitCollection> hf;
0153 iEvent.getByToken(hfToken_, hf);
0154 if (!hf.isValid()) {
0155
0156 if (!allowMissingInputs_) {
0157 *hf;
0158 }
0159 } else {
0160 for (auto const hfItr : (HFRecHitCollection)(*hf)) {
0161 DetId id = hfItr.detid();
0162 theFillDetMap0_[id] += 1.;
0163 theFillDetMap1_[id] += hfItr.energy();
0164 theFillDetMap2_[id] += pow(hfItr.energy(), 2);
0165 theFillDetMap3_[id] += pow(hfItr.energy(), 3);
0166 theFillDetMap4_[id] += pow(hfItr.energy(), 4);
0167 }
0168 }
0169 }
0170 }
0171 }