RecoJets/JetAnalyzers/src/JetToDigiDump.cc - master - githubmirror/cmssw

RecoJets/JetAnalyzers/src/JetToDigiDump.cc
349 строк20 KB
Cms Build
Clang-Format
03 июн 2019, 06:54
03 июн 2019, 06:54c1e01f0
// JetToDigiDump.cc
// Description:  Prints out Jets, consituent CaloTowers, constituent RecHits and associated Digis (the digis for HCAL only).
//               The user can specify which level in the config file:
//               DumpLevel="Jets":    Printout of jets and their kinematic quantities.
//               DumpLevel="Towers":  Nested Printout of jets and their constituent CaloTowers
//               DumpLevel="RecHits": Nested Printout of jets, constituent CaloTowers and constituent RecHits
//               DumpLevel="Digis":   Nested Printout of jets, constituent CaloTowers, RecHits and all the HCAL digis
//                                    associated with the RecHit channel (no links exist to go back to actual digis used).
//               Does simple sanity checks on energy sums at each level: jets=sum of towers, tower=sum of RecHits.
//               Does quick and dirty estimate of the fC/GeV factor that was applied to make the RecHit from the Digis.
//
// Author: Robert M. Harris
// Date:  19 - October - 2006
//
#include "RecoJets/JetAnalyzers/interface/JetToDigiDump.h"
#include "DataFormats/JetReco/interface/CaloJetCollection.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
//in CaloJet: #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
// just for the CaloTowerPtr declaration:
// in CaloTowerCollection: #include "DataFormats/CaloTowers/interface/CaloTower.h"
// in CaloTowerCollection: #include "DataFormats/CaloTowers/interface/CaloTowerDefs.h"
#include "DataFormats/DetId/interface/DetId.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
#include "DataFormats/EcalDetId/interface/EBDetId.h"
#include "DataFormats/EcalDetId/interface/EEDetId.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
#include "DataFormats/EcalDigi/interface/EBDataFrame.h"
#include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include <TROOT.h>
#include <TSystem.h>
#include <TFile.h>
#include <TCanvas.h>
#include <cmath>
using namespace edm;
using namespace reco;
using namespace std;

JetToDigiDump::JetToDigiDump(const ParameterSet& cfg)
    : DumpLevel(cfg.getParameter<string>("DumpLevel")),
      CaloJetAlg(cfg.getParameter<string>("CaloJetAlg")),
      DebugLevel(cfg.getParameter<int>("DebugLevel")),
      ShowECal(cfg.getParameter<bool>("ShowECal")) {}

void JetToDigiDump::beginJob() {
  if (DumpLevel == "Jets") {
    cout << "Dump of Jets" << endl;
    Dump = 1;
  } else if (DumpLevel == "Towers") {
    cout << "Dump of Jets and constituent CaloTowers" << endl;
    Dump = 2;
  } else if (DumpLevel == "RecHits") {
    cout << "Dump of Jets, constituent CaloTowers, and constituent RecHits" << endl;
    Dump = 3;
  } else if (DumpLevel == "Digis") {
    cout << "Dump of Jets, constituent CaloTowers, constituent RecHits and associated Digis" << endl;
    Dump = 4;
  }
  cout << "Jet Algorithm being dumped is " << CaloJetAlg << endl;
  cout << "Debug level is " << DebugLevel << endl;
  //Initialize some stuff
  evtCount = 0;
}

void JetToDigiDump::analyze(const Event& evt, const EventSetup& es) {
  int jetInd;
  Handle<CaloJetCollection> caloJets;
  Handle<CaloTowerCollection> caloTowers;
  Handle<HBHERecHitCollection> HBHERecHits;
  Handle<HORecHitCollection> HORecHits;
  Handle<HFRecHitCollection> HFRecHits;
  Handle<EBRecHitCollection> EBRecHits;
  Handle<EERecHitCollection> EERecHits;
  Handle<HBHEDigiCollection> HBHEDigis;
  Handle<HODigiCollection> HODigis;
  Handle<HFDigiCollection> HFDigis;
  Handle<EEDigiCollection> EEDigis;
  Handle<EBDigiCollection> EBDigis;
  // Old:
  //Handle<edm::SortedCollection<EBDataFrame> > EBDigis;
  //Handle<edm::SortedCollection<EBDataFrame> > EEDigis;

  //Find the CaloTowers in leading CaloJets
  if (DebugLevel)
    cout << "Getting caloJets" << endl;

  evt.getByLabel(CaloJetAlg, caloJets);
  if (Dump >= 2)
    evt.getByLabel("towerMaker", caloTowers);
  if (Dump >= 3) {
    if (DebugLevel)
      cout << "Getting recHits" << endl;
    evt.getByLabel("hbhereco", HBHERecHits);
    evt.getByLabel("horeco", HORecHits);
    evt.getByLabel("hfreco", HFRecHits);
    evt.getByLabel("ecalRecHit", "EcalRecHitsEB", EBRecHits);
    evt.getByLabel("ecalRecHit", "EcalRecHitsEE", EERecHits);
    if (DebugLevel)
      cout << "# of hits gotten - HBHE: " << HBHERecHits->size() << endl;
    evt.getByLabel("hcalDigis", HBHEDigis);
    evt.getByLabel("hcalDigis", HODigis);
    evt.getByLabel("hcalDigis", HFDigis);
    if (DebugLevel)
      cout << "# of digis gotten - HBHE: " << HBHEDigis->size() << endl;
    if (ShowECal) {
      evt.getByLabel("ecalDigis", "ebDigis", EBDigis);
      evt.getByLabel("ecalDigis", "eeDigis", EEDigis);
    }
  }

  cout << endl << "Evt: " << evtCount << ", Num Jets=" << caloJets->end() - caloJets->begin() << endl;
  if (Dump >= 1)
    cout << "   *********************************************************" << endl;
  jetInd = 0;
  if (Dump >= 1)
    for (CaloJetCollection::const_iterator jet = caloJets->begin(); jet != caloJets->end(); ++jet) {
      //2_1_?    std::vector <CaloTowerPtr> towers = jet->getCaloConstituents ();
      //2_0_7"
      std::vector<CaloTowerPtr> towers = jet->getCaloConstituents();
      int nConstituents = towers.size();
      cout << "   Jet: " << jetInd << ", eta=" << jet->eta() << ", phi=" << jet->phi() << ", pt=" << jet->pt()
           << ",E=" << jet->energy() << ", EB E=" << jet->emEnergyInEB() << " ,HB E=" << jet->hadEnergyInHB()
           << ", HO E=" << jet->hadEnergyInHO() << " ,EE E=" << jet->emEnergyInEE() << ", HE E=" << jet->hadEnergyInHE()
           << ", HF E=" << jet->hadEnergyInHF() + jet->emEnergyInHF() << ", Num Towers=" << nConstituents << endl;
      if (Dump >= 2)
        cout << "      =====================================================" << endl;
      float sumTowerE = 0.0;
      if (Dump >= 2)
        for (int i = 0; i < nConstituents; i++) {
          CaloTowerCollection::const_iterator theTower =
              caloTowers->find(towers[i]->id());  //Find the tower from its CaloTowerDetID
          if (theTower == caloTowers->end()) {
            cerr << "Bug? Can't find the tower" << endl;
            return;
          }
          int ietaTower = towers[i]->id().ieta();
          int iphiTower = towers[i]->id().iphi();
          sumTowerE += theTower->energy();
          size_t numRecHits = theTower->constituentsSize();
          cout << "      Tower " << i << ": ieta=" << ietaTower << ", eta=" << theTower->eta() << ", iphi=" << iphiTower
               << ", phi=" << theTower->phi() << ", energy=" << theTower->energy() << ", EM=" << theTower->emEnergy()
               << ", HAD=" << theTower->hadEnergy() << ", HO=" << theTower->outerEnergy()
               << ",  Num Rec Hits =" << numRecHits << endl;
          if (Dump >= 3)
            cout << "         ------------------------------------------------" << endl;
          float sumRecHitE = 0.0;
          if (Dump >= 3)
            for (size_t j = 0; j < numRecHits; j++) {
              DetId RecHitDetID = theTower->constituent(j);
              DetId::Detector DetNum = RecHitDetID.det();
              if (DetNum == DetId::Hcal) {
                //cout << "RecHit " << j << ": Detector = " << DetNum << ": Hcal " << endl;
                HcalDetId HcalID = RecHitDetID;
                HcalSubdetector HcalNum = HcalID.subdet();
                if (HcalNum == HcalBarrel) {
                  HBHERecHitCollection::const_iterator theRecHit = HBHERecHits->find(HcalID);
                  sumRecHitE += theRecHit->energy();
                  HBHEDigiCollection::const_iterator theDigis = HBHEDigis->find(HcalID);
                  cout << "         RecHit: " << j << ": HB, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
                       << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
                       << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
                       << ", presamples =" << theDigis->presamples() << endl;
                  /* const HcalElectronicsId HW_ID = theDigis->elecId();
                cout << "Digi: Index=" << HW_ID.linearIndex() << ", raw ID=" <<  HW_ID.rawId() << ", fiberChan=" << HW_ID.fiberChanId() <<  ", fiberInd=" <<  HW_ID.fiberIndex() \
	             << ", HTR chan=" <<  HW_ID.htrChanId() << ", HTR Slot=" <<   HW_ID.htrSlot() << ", HDR top/bot=" << HW_ID.htrTopBottom() \
	        << ", VME crate=" <<  HW_ID.readoutVMECrateId() << ", DCC=" << HW_ID.dccid() << ", DCC spigot=" <<  HW_ID.spigot() << endl;
              */
                  float SumDigiCharge = 0.0;
                  float EstimatedPedestal = 0.0;
                  int SamplesToAdd = 4;
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                  if (Dump >= 4)
                    for (int k = 0; k < theDigis->size(); k++) {
                      const HcalQIESample QIE = theDigis->sample(k);
                      if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
                        SumDigiCharge += QIE.nominal_fC();
                      if (k < theDigis->presamples() - 1)
                        EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
                      cout << "            Digi: " << k << ", cap ID = " << QIE.capid()
                           << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
                    }
                  if (Dump >= 4)
                    cout << "            4 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
                         << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                } else if (HcalNum == HcalEndcap) {
                  HBHERecHitCollection::const_iterator theRecHit = HBHERecHits->find(HcalID);
                  if ((abs(HcalID.ieta()) == 28 || abs(HcalID.ieta()) == 29) && HcalID.depth() == 3) {
                    sumRecHitE += theRecHit->energy() / 2;  //Depth 3 split over tower 28 & 29
                  } else {
                    sumRecHitE += theRecHit->energy();
                  }
                  HBHEDigiCollection::const_iterator theDigis = HBHEDigis->find(HcalID);
                  cout << "         RecHit: " << j << ": HE, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
                       << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
                       << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
                       << ", presamples =" << theDigis->presamples() << endl;
                  float SumDigiCharge = 0.0;
                  float EstimatedPedestal = 0.0;
                  int SamplesToAdd = 4;
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                  if (Dump >= 4)
                    for (int k = 0; k < theDigis->size(); k++) {
                      const HcalQIESample QIE = theDigis->sample(k);
                      if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
                        SumDigiCharge += QIE.nominal_fC();
                      if (k < theDigis->presamples() - 1)
                        EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
                      cout << "            Digi: " << k << ", cap ID = " << QIE.capid()
                           << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
                    }
                  if (Dump >= 4)
                    cout << "            4 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
                         << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                } else if (HcalNum == HcalOuter) {
                  HORecHitCollection::const_iterator theRecHit = HORecHits->find(HcalID);
                  sumRecHitE += theRecHit->energy();
                  HODigiCollection::const_iterator theDigis = HODigis->find(HcalID);
                  cout << "         RecHit: " << j << ": HO, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
                       << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
                       << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
                       << ", presamples =" << theDigis->presamples() << endl;
                  float SumDigiCharge = 0.0;
                  float EstimatedPedestal = 0.0;
                  int SamplesToAdd = 4;
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                  if (Dump >= 4)
                    for (int k = 0; k < theDigis->size(); k++) {
                      const HcalQIESample QIE = theDigis->sample(k);
                      if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
                        SumDigiCharge += QIE.nominal_fC();
                      if (k < theDigis->presamples() - 1)
                        EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
                      cout << "            Digi: " << k << ", cap ID = " << QIE.capid()
                           << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
                    }
                  if (Dump >= 4)
                    cout << "            4 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
                         << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                } else if (HcalNum == HcalForward) {
                  HFRecHitCollection::const_iterator theRecHit = HFRecHits->find(HcalID);
                  sumRecHitE += theRecHit->energy();
                  HFDigiCollection::const_iterator theDigis = HFDigis->find(HcalID);
                  cout << "         RecHit: " << j << ": HF, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
                       << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
                       << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
                       << ", presamples =" << theDigis->presamples() << endl;
                  float SumDigiCharge = 0.0;
                  float EstimatedPedestal = 0.0;
                  int SamplesToAdd = 1;
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                  if (Dump >= 4)
                    for (int k = 0; k < theDigis->size(); k++) {
                      const HcalQIESample QIE = theDigis->sample(k);
                      if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
                        SumDigiCharge += QIE.nominal_fC();
                      if (k < theDigis->presamples() - 1)
                        EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
                      cout << "            Digi: " << k << ", cap ID = " << QIE.capid()
                           << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
                    }
                  if (Dump >= 4)
                    cout << "            1 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
                         << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                }
              }
              if (ShowECal && DetNum == DetId::Ecal) {
                int EcalNum = RecHitDetID.subdetId();
                if (EcalNum == 1) {
                  EBDetId EcalID = RecHitDetID;
                  EBRecHitCollection::const_iterator theRecHit = EBRecHits->find(EcalID);
                  EBDigiCollection::const_iterator theDigis = EBDigis->find(EcalID);
                  sumRecHitE += theRecHit->energy();
                  cout << "         RecHit " << j << ": EB, ieta=" << EcalID.ieta() << ", iphi=" << EcalID.iphi()
                       << ", SM=" << EcalID.ism() << ", energy=" << theRecHit->energy()
                       << ", All Digis=" << theDigis->size() << endl;
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                  if (Dump >= 4)
                    for (unsigned int k = 0; k < theDigis->size(); k++) {
                      EBDataFrame frame(*theDigis);
                      const EcalMGPASample MGPA = frame.sample(k);
                      cout << "            Digi: " << k << ": ADC Sample = " << MGPA.adc()
                           << ", Gain ID = " << MGPA.gainId() << endl;
                    }
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                } else if (EcalNum == 2) {
                  EEDetId EcalID = RecHitDetID;
                  EERecHitCollection::const_iterator theRecHit = EERecHits->find(EcalID);
                  EEDigiCollection::const_iterator theDigis = EEDigis->find(EcalID);
                  sumRecHitE += theRecHit->energy();
                  cout << "         RecHit " << j << ": EE, ix=" << EcalID.ix() << ", iy=" << EcalID.iy()
                       << ", energy=" << theRecHit->energy() << ", All Digis=" << theDigis->size() << endl;
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                  if (Dump >= 4)
                    for (unsigned int k = 0; k < theDigis->size(); k++) {
                      EEDataFrame frame(*theDigis);
                      const EcalMGPASample MGPA = frame.sample(k);
                      cout << "            Digi: " << k << ": ADC Sample = " << MGPA.adc()
                           << ", Gain ID = " << MGPA.gainId() << endl;
                    }
                  if (Dump >= 4)
                    cout << "            ......................................" << endl;
                }
              }
            }
          if (Dump >= 3) {
            if (abs(ietaTower) == 28 || abs(ietaTower) == 29) {
              cout << "         Splitted Sum of RecHit Energies=" << sumRecHitE
                   << ", CaloTower energy=" << theTower->energy() << endl;
            } else {
              cout << "         Sum of RecHit Energies=" << sumRecHitE << ", CaloTower energy=" << theTower->energy()
                   << endl;
            }
          }
          if (Dump >= 3)
            cout << "         ------------------------------------------------" << endl;
        }
      if (Dump >= 2)
        cout << "      Sum of tower energies=" << sumTowerE << ", CaloJet energy=" << jet->energy() << endl;
      jetInd++;
      if (Dump >= 2)
        cout << "      =====================================================" << endl;
    }
  evtCount++;
  if (Dump >= 1)
    cout << "   *********************************************************" << endl;
}

void JetToDigiDump::endJob() {}
#include "FWCore/Framework/interface/MakerMacros.h"
DEFINE_FWK_MODULE(JetToDigiDump);