MinBias.cc

Go to the documentation of this file.

// user include files
#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "JetMETCorrections/MinBias/interface/MinBias.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "Geometry/HcalTowerAlgo/interface/HcalGeometry.h"

namespace cms {
  MinBias::MinBias(const edm::ParameterSet& iConfig) : geo_(nullptr) {
    // get names of modules, producing object collections
    hbheLabel_ = iConfig.getParameter<std::string>("hbheInput");
    hoLabel_ = iConfig.getParameter<std::string>("hoInput");
    hfLabel_ = iConfig.getParameter<std::string>("hfInput");
    hbheToken_ = mayConsume<HBHERecHitCollection>(edm::InputTag(hbheLabel_));
    hoToken_ = mayConsume<HORecHitCollection>(edm::InputTag(hoLabel_));
    hfToken_ = mayConsume<HFRecHitCollection>(edm::InputTag(hfLabel_));
    allowMissingInputs_ = iConfig.getUntrackedParameter<bool>("AllowMissingInputs", false);
  }

  void MinBias::beginJob() {
    edm::Service<TFileService> fs;
    myTree_ = fs->make<TTree>("RecJet", "RecJet Tree");
    myTree_->Branch("mydet", &mydet, "mydet/I");
    myTree_->Branch("mysubd", &mysubd, "mysubd/I");
    myTree_->Branch("depth", &depth, "depth/I");
    myTree_->Branch("ieta", &ieta, "ieta/I");
    myTree_->Branch("iphi", &iphi, "iphi/I");
    myTree_->Branch("eta", &eta, "eta/F");
    myTree_->Branch("phi", &phi, "phi/F");
    myTree_->Branch("mom1", &mom1, "mom1/F");
    myTree_->Branch("mom2", &mom2, "mom2/F");
    myTree_->Branch("mom3", &mom3, "mom3/F");
    myTree_->Branch("mom4", &mom4, "mom4/F");
  }

  void MinBias::beginRun(edm::Run const& iRun, edm::EventSetup const& iSetup) {
    edm::ESHandle<CaloGeometry> pG;
    iSetup.get<CaloGeometryRecord>().get(pG);
    geo_ = pG.product();
    std::vector<DetId> did = geo_->getValidDetIds();

    for (auto const& id : did) {
      if ((id).det() == DetId::Hcal) {
        theFillDetMap0_[id] = 0.;
        theFillDetMap1_[id] = 0.;
        theFillDetMap2_[id] = 0.;
        theFillDetMap3_[id] = 0.;
        theFillDetMap4_[id] = 0.;
      }
    }
  }

  void MinBias::endJob() {
    const HcalGeometry* hgeo = static_cast<const HcalGeometry*>(geo_->getSubdetectorGeometry(DetId::Hcal, 1));
    const std::vector<DetId>& did = hgeo->getValidDetIds();
    int i = 0;
    for (const auto& id : did) {
      //      if( id.det() == DetId::Hcal ) {
      GlobalPoint pos = hgeo->getPosition(id);
      mydet = (int)(id.det());
      mysubd = (id.subdetId());
      depth = HcalDetId(id).depth();
      ieta = HcalDetId(id).ieta();
      iphi = HcalDetId(id).iphi();
      phi = pos.phi();
      eta = pos.eta();
      if (theFillDetMap0_[id] > 0.) {
        mom1 = theFillDetMap1_[id] / theFillDetMap0_[id];
        mom2 = theFillDetMap2_[id] / theFillDetMap0_[id] - (mom1 * mom1);
        mom3 = theFillDetMap3_[id] / theFillDetMap0_[id] - 3. * mom1 * theFillDetMap2_[id] / theFillDetMap0_[id] +
               2. * pow(mom2, 3);
        mom4 = (theFillDetMap4_[id] - 4. * mom1 * theFillDetMap3_[id] + 6. * pow(mom1, 2) * theFillDetMap2_[id]) /
                   theFillDetMap0_[id] -
               3. * pow(mom1, 4);

      } else {
        mom1 = 0.;
        mom2 = 0.;
        mom3 = 0.;
        mom4 = 0.;
      }
      edm::LogWarning("MinBias") << " Detector " << id.rawId() << " mydet " << mydet << " " << mysubd << " " << depth
                                 << " " << ieta << " " << iphi << " " << pos.eta() << " " << pos.phi();
      edm::LogWarning("MinBias") << " Energy " << mom1 << " " << mom2 << std::endl;
      myTree_->Fill();
      i++;
      //      }
    }
    edm::LogWarning("MinBias") << " The number of CaloDet records " << did.size();
    edm::LogWarning("MinBias") << " The number of Hcal records " << i << std::endl;

    /*
  std::cout << "===== Start writing user histograms =====" << std::endl;
  hOutputFile->SetCompressionLevel(2);
  hOutputFile->cd();
  myTree->Write();
  hOutputFile->Close() ;
  std::cout << "===== End writing user histograms =======" << std::endl;
  */
  }

  //
  // member functions
  //

  // ------------ method called to produce the data  ------------
  void MinBias::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
    if (!hbheLabel_.empty()) {
      edm::Handle<HBHERecHitCollection> hbhe;
      iEvent.getByToken(hbheToken_, hbhe);
      if (!hbhe.isValid()) {
        // can't find it!
        if (!allowMissingInputs_) {
          *hbhe;  // will throw the proper exception
        }
      } else {
        for (auto const& hbheItr : (HBHERecHitCollection)(*hbhe)) {
          DetId id = (hbheItr).detid();
          if (hbheItr.energy() > 0.)
            edm::LogWarning("MinBias") << " Energy = " << hbheItr.energy();
          theFillDetMap0_[id] += 1.;
          theFillDetMap1_[id] += hbheItr.energy();
          theFillDetMap2_[id] += pow(hbheItr.energy(), 2);
          theFillDetMap3_[id] += pow(hbheItr.energy(), 3);
          theFillDetMap4_[id] += pow(hbheItr.energy(), 4);
        }
      }
    }

    if (!hoLabel_.empty()) {
      edm::Handle<HORecHitCollection> ho;
      iEvent.getByToken(hoToken_, ho);
      if (!ho.isValid()) {
        // can't find it!
        if (!allowMissingInputs_) {
          *ho;  // will throw the proper exception
        }
      } else {
        for (auto const& hoItr : (HORecHitCollection)(*ho)) {
          DetId id = hoItr.detid();
          theFillDetMap0_[id] += 1.;
          theFillDetMap1_[id] += hoItr.energy();
          theFillDetMap2_[id] += pow(hoItr.energy(), 2);
          theFillDetMap3_[id] += pow(hoItr.energy(), 3);
          theFillDetMap4_[id] += pow(hoItr.energy(), 4);
        }
      }
    }

    if (!hfLabel_.empty()) {
      edm::Handle<HFRecHitCollection> hf;
      iEvent.getByToken(hfToken_, hf);
      if (!hf.isValid()) {
        // can't find it!
        if (!allowMissingInputs_) {
          *hf;  // will throw the proper exception
        }
      } else {
        for (auto const hfItr : (HFRecHitCollection)(*hf)) {
          DetId id = hfItr.detid();
          theFillDetMap0_[id] += 1.;
          theFillDetMap1_[id] += hfItr.energy();
          theFillDetMap2_[id] += pow(hfItr.energy(), 2);
          theFillDetMap3_[id] += pow(hfItr.energy(), 3);
          theFillDetMap4_[id] += pow(hfItr.energy(), 4);
        }
      }
    }
  }
}  // namespace cms