MinBias/src/MinBias.cc

0001 // user include files
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     // get names of modules, producing object collections
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       //      if( id.det() == DetId::Hcal ) {
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   std::cout << "===== Start writing user histograms =====" << std::endl;
0095   hOutputFile->SetCompressionLevel(2);
0096   hOutputFile->cd();
0097   myTree->Write();
0098   hOutputFile->Close() ;
0099   std::cout << "===== End writing user histograms =======" << std::endl;
0100   */
0101   }
0102
0103   //
0104   // member functions
0105   //
0106
0107   // ------------ method called to produce the data  ------------
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         // can't find it!
0114         if (!allowMissingInputs_) {
0115           *hbhe;  // will throw the proper exception
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         // can't find it!
0136         if (!allowMissingInputs_) {
0137           *ho;  // will throw the proper exception
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         // can't find it!
0156         if (!allowMissingInputs_) {
0157           *hf;  // will throw the proper exception
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 }  // namespace cms