Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:27 

0001 // JetToDigiDump.cc

0002 // Description:  Prints out Jets, consituent CaloTowers, constituent RecHits and associated Digis (the digis for HCAL only).

0003 //               The user can specify which level in the config file:

0004 //               DumpLevel="Jets":    Printout of jets and their kinematic quantities.

0005 //               DumpLevel="Towers":  Nested Printout of jets and their constituent CaloTowers

0006 //               DumpLevel="RecHits": Nested Printout of jets, constituent CaloTowers and constituent RecHits

0007 //               DumpLevel="Digis":   Nested Printout of jets, constituent CaloTowers, RecHits and all the HCAL digis

0008 //                                    associated with the RecHit channel (no links exist to go back to actual digis used).

0009 //               Does simple sanity checks on energy sums at each level: jets=sum of towers, tower=sum of RecHits.

0010 //               Does quick and dirty estimate of the fC/GeV factor that was applied to make the RecHit from the Digis.

0011 //

0012 // Author: Robert M. Harris

0013 // Date:  19 - October - 2006

0014 //

0015 #include "RecoJets/JetAnalyzers/interface/JetToDigiDump.h"
0016 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
0017 #include "DataFormats/JetReco/interface/CaloJet.h"
0018 //in CaloJet: #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"

0019 // just for the CaloTowerPtr declaration:

0020 // in CaloTowerCollection: #include "DataFormats/CaloTowers/interface/CaloTower.h"

0021 // in CaloTowerCollection: #include "DataFormats/CaloTowers/interface/CaloTowerDefs.h"

0022 #include "DataFormats/DetId/interface/DetId.h"
0023 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
0024 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
0025 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
0026 #include "DataFormats/EcalDetId/interface/EBDetId.h"
0027 #include "DataFormats/EcalDetId/interface/EEDetId.h"
0028 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
0029 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
0030 #include "DataFormats/EcalDigi/interface/EBDataFrame.h"
0031 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
0032 #include "FWCore/Framework/interface/Event.h"
0033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0034 #include <TROOT.h>
0035 #include <TSystem.h>
0036 #include <TFile.h>
0037 #include <TCanvas.h>
0038 #include <cmath>
0039 using namespace edm;
0040 using namespace reco;
0041 using namespace std;
0042 
0043 JetToDigiDump::JetToDigiDump(const ParameterSet& cfg)
0044     : DumpLevel(cfg.getParameter<string>("DumpLevel")),
0045       CaloJetAlg(cfg.getParameter<string>("CaloJetAlg")),
0046       DebugLevel(cfg.getParameter<int>("DebugLevel")),
0047       ShowECal(cfg.getParameter<bool>("ShowECal")) {}
0048 
0049 void JetToDigiDump::beginJob() {
0050   if (DumpLevel == "Jets") {
0051     cout << "Dump of Jets" << endl;
0052     Dump = 1;
0053   } else if (DumpLevel == "Towers") {
0054     cout << "Dump of Jets and constituent CaloTowers" << endl;
0055     Dump = 2;
0056   } else if (DumpLevel == "RecHits") {
0057     cout << "Dump of Jets, constituent CaloTowers, and constituent RecHits" << endl;
0058     Dump = 3;
0059   } else if (DumpLevel == "Digis") {
0060     cout << "Dump of Jets, constituent CaloTowers, constituent RecHits and associated Digis" << endl;
0061     Dump = 4;
0062   }
0063   cout << "Jet Algorithm being dumped is " << CaloJetAlg << endl;
0064   cout << "Debug level is " << DebugLevel << endl;
0065   //Initialize some stuff

0066   evtCount = 0;
0067 }
0068 
0069 void JetToDigiDump::analyze(const Event& evt, const EventSetup& es) {
0070   int jetInd;
0071   Handle<CaloJetCollection> caloJets;
0072   Handle<CaloTowerCollection> caloTowers;
0073   Handle<HBHERecHitCollection> HBHERecHits;
0074   Handle<HORecHitCollection> HORecHits;
0075   Handle<HFRecHitCollection> HFRecHits;
0076   Handle<EBRecHitCollection> EBRecHits;
0077   Handle<EERecHitCollection> EERecHits;
0078   Handle<HBHEDigiCollection> HBHEDigis;
0079   Handle<HODigiCollection> HODigis;
0080   Handle<HFDigiCollection> HFDigis;
0081   Handle<EEDigiCollection> EEDigis;
0082   Handle<EBDigiCollection> EBDigis;
0083   // Old:

0084   //Handle<edm::SortedCollection<EBDataFrame> > EBDigis;

0085   //Handle<edm::SortedCollection<EBDataFrame> > EEDigis;

0086 
0087   //Find the CaloTowers in leading CaloJets

0088   if (DebugLevel)
0089     cout << "Getting caloJets" << endl;
0090 
0091   evt.getByLabel(CaloJetAlg, caloJets);
0092   if (Dump >= 2)
0093     evt.getByLabel("towerMaker", caloTowers);
0094   if (Dump >= 3) {
0095     if (DebugLevel)
0096       cout << "Getting recHits" << endl;
0097     evt.getByLabel("hbhereco", HBHERecHits);
0098     evt.getByLabel("horeco", HORecHits);
0099     evt.getByLabel("hfreco", HFRecHits);
0100     evt.getByLabel("ecalRecHit", "EcalRecHitsEB", EBRecHits);
0101     evt.getByLabel("ecalRecHit", "EcalRecHitsEE", EERecHits);
0102     if (DebugLevel)
0103       cout << "# of hits gotten - HBHE: " << HBHERecHits->size() << endl;
0104     evt.getByLabel("hcalDigis", HBHEDigis);
0105     evt.getByLabel("hcalDigis", HODigis);
0106     evt.getByLabel("hcalDigis", HFDigis);
0107     if (DebugLevel)
0108       cout << "# of digis gotten - HBHE: " << HBHEDigis->size() << endl;
0109     if (ShowECal) {
0110       evt.getByLabel("ecalDigis", "ebDigis", EBDigis);
0111       evt.getByLabel("ecalDigis", "eeDigis", EEDigis);
0112     }
0113   }
0114 
0115   cout << endl << "Evt: " << evtCount << ", Num Jets=" << caloJets->end() - caloJets->begin() << endl;
0116   if (Dump >= 1)
0117     cout << "   *********************************************************" << endl;
0118   jetInd = 0;
0119   if (Dump >= 1)
0120     for (CaloJetCollection::const_iterator jet = caloJets->begin(); jet != caloJets->end(); ++jet) {
0121       //2_1_?    std::vector <CaloTowerPtr> towers = jet->getCaloConstituents ();

0122       //2_0_7"

0123       std::vector<CaloTowerPtr> towers = jet->getCaloConstituents();
0124       int nConstituents = towers.size();
0125       cout << "   Jet: " << jetInd << ", eta=" << jet->eta() << ", phi=" << jet->phi() << ", pt=" << jet->pt()
0126            << ",E=" << jet->energy() << ", EB E=" << jet->emEnergyInEB() << " ,HB E=" << jet->hadEnergyInHB()
0127            << ", HO E=" << jet->hadEnergyInHO() << " ,EE E=" << jet->emEnergyInEE() << ", HE E=" << jet->hadEnergyInHE()
0128            << ", HF E=" << jet->hadEnergyInHF() + jet->emEnergyInHF() << ", Num Towers=" << nConstituents << endl;
0129       if (Dump >= 2)
0130         cout << "      =====================================================" << endl;
0131       float sumTowerE = 0.0;
0132       if (Dump >= 2)
0133         for (int i = 0; i < nConstituents; i++) {
0134           CaloTowerCollection::const_iterator theTower =
0135               caloTowers->find(towers[i]->id());  //Find the tower from its CaloTowerDetID

0136           if (theTower == caloTowers->end()) {
0137             cerr << "Bug? Can't find the tower" << endl;
0138             return;
0139           }
0140           int ietaTower = towers[i]->id().ieta();
0141           int iphiTower = towers[i]->id().iphi();
0142           sumTowerE += theTower->energy();
0143           size_t numRecHits = theTower->constituentsSize();
0144           cout << "      Tower " << i << ": ieta=" << ietaTower << ", eta=" << theTower->eta() << ", iphi=" << iphiTower
0145                << ", phi=" << theTower->phi() << ", energy=" << theTower->energy() << ", EM=" << theTower->emEnergy()
0146                << ", HAD=" << theTower->hadEnergy() << ", HO=" << theTower->outerEnergy()
0147                << ",  Num Rec Hits =" << numRecHits << endl;
0148           if (Dump >= 3)
0149             cout << "         ------------------------------------------------" << endl;
0150           float sumRecHitE = 0.0;
0151           if (Dump >= 3)
0152             for (size_t j = 0; j < numRecHits; j++) {
0153               DetId RecHitDetID = theTower->constituent(j);
0154               DetId::Detector DetNum = RecHitDetID.det();
0155               if (DetNum == DetId::Hcal) {
0156                 //cout << "RecHit " << j << ": Detector = " << DetNum << ": Hcal " << endl;

0157                 HcalDetId HcalID = RecHitDetID;
0158                 HcalSubdetector HcalNum = HcalID.subdet();
0159                 if (HcalNum == HcalBarrel) {
0160                   HBHERecHitCollection::const_iterator theRecHit = HBHERecHits->find(HcalID);
0161                   sumRecHitE += theRecHit->energy();
0162                   HBHEDigiCollection::const_iterator theDigis = HBHEDigis->find(HcalID);
0163                   cout << "         RecHit: " << j << ": HB, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
0164                        << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
0165                        << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
0166                        << ", presamples =" << theDigis->presamples() << endl;
0167                   /* const HcalElectronicsId HW_ID = theDigis->elecId();

0168                 cout << "Digi: Index=" << HW_ID.linearIndex() << ", raw ID=" <<  HW_ID.rawId() << ", fiberChan=" << HW_ID.fiberChanId() <<  ", fiberInd=" <<  HW_ID.fiberIndex() \

0169                  << ", HTR chan=" <<  HW_ID.htrChanId() << ", HTR Slot=" <<   HW_ID.htrSlot() << ", HDR top/bot=" << HW_ID.htrTopBottom() \

0170             << ", VME crate=" <<  HW_ID.readoutVMECrateId() << ", DCC=" << HW_ID.dccid() << ", DCC spigot=" <<  HW_ID.spigot() << endl;

0171               */
0172                   float SumDigiCharge = 0.0;
0173                   float EstimatedPedestal = 0.0;
0174                   int SamplesToAdd = 4;
0175                   if (Dump >= 4)
0176                     cout << "            ......................................" << endl;
0177                   if (Dump >= 4)
0178                     for (int k = 0; k < theDigis->size(); k++) {
0179                       const HcalQIESample QIE = theDigis->sample(k);
0180                       if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
0181                         SumDigiCharge += QIE.nominal_fC();
0182                       if (k < theDigis->presamples() - 1)
0183                         EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
0184                       cout << "            Digi: " << k << ", cap ID = " << QIE.capid()
0185                            << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
0186                     }
0187                   if (Dump >= 4)
0188                     cout << "            4 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
0189                          << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
0190                   if (Dump >= 4)
0191                     cout << "            ......................................" << endl;
0192                 } else if (HcalNum == HcalEndcap) {
0193                   HBHERecHitCollection::const_iterator theRecHit = HBHERecHits->find(HcalID);
0194                   if ((abs(HcalID.ieta()) == 28 || abs(HcalID.ieta()) == 29) && HcalID.depth() == 3) {
0195                     sumRecHitE += theRecHit->energy() / 2;  //Depth 3 split over tower 28 & 29

0196                   } else {
0197                     sumRecHitE += theRecHit->energy();
0198                   }
0199                   HBHEDigiCollection::const_iterator theDigis = HBHEDigis->find(HcalID);
0200                   cout << "         RecHit: " << j << ": HE, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
0201                        << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
0202                        << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
0203                        << ", presamples =" << theDigis->presamples() << endl;
0204                   float SumDigiCharge = 0.0;
0205                   float EstimatedPedestal = 0.0;
0206                   int SamplesToAdd = 4;
0207                   if (Dump >= 4)
0208                     cout << "            ......................................" << endl;
0209                   if (Dump >= 4)
0210                     for (int k = 0; k < theDigis->size(); k++) {
0211                       const HcalQIESample QIE = theDigis->sample(k);
0212                       if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
0213                         SumDigiCharge += QIE.nominal_fC();
0214                       if (k < theDigis->presamples() - 1)
0215                         EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
0216                       cout << "            Digi: " << k << ", cap ID = " << QIE.capid()
0217                            << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
0218                     }
0219                   if (Dump >= 4)
0220                     cout << "            4 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
0221                          << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
0222                   if (Dump >= 4)
0223                     cout << "            ......................................" << endl;
0224                 } else if (HcalNum == HcalOuter) {
0225                   HORecHitCollection::const_iterator theRecHit = HORecHits->find(HcalID);
0226                   sumRecHitE += theRecHit->energy();
0227                   HODigiCollection::const_iterator theDigis = HODigis->find(HcalID);
0228                   cout << "         RecHit: " << j << ": HO, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
0229                        << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
0230                        << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
0231                        << ", presamples =" << theDigis->presamples() << endl;
0232                   float SumDigiCharge = 0.0;
0233                   float EstimatedPedestal = 0.0;
0234                   int SamplesToAdd = 4;
0235                   if (Dump >= 4)
0236                     cout << "            ......................................" << endl;
0237                   if (Dump >= 4)
0238                     for (int k = 0; k < theDigis->size(); k++) {
0239                       const HcalQIESample QIE = theDigis->sample(k);
0240                       if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
0241                         SumDigiCharge += QIE.nominal_fC();
0242                       if (k < theDigis->presamples() - 1)
0243                         EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
0244                       cout << "            Digi: " << k << ", cap ID = " << QIE.capid()
0245                            << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
0246                     }
0247                   if (Dump >= 4)
0248                     cout << "            4 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
0249                          << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
0250                   if (Dump >= 4)
0251                     cout << "            ......................................" << endl;
0252                 } else if (HcalNum == HcalForward) {
0253                   HFRecHitCollection::const_iterator theRecHit = HFRecHits->find(HcalID);
0254                   sumRecHitE += theRecHit->energy();
0255                   HFDigiCollection::const_iterator theDigis = HFDigis->find(HcalID);
0256                   cout << "         RecHit: " << j << ": HF, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()
0257                        << ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy()
0258                        << ", time=" << theRecHit->time() << ", All Digis=" << theDigis->size()
0259                        << ", presamples =" << theDigis->presamples() << endl;
0260                   float SumDigiCharge = 0.0;
0261                   float EstimatedPedestal = 0.0;
0262                   int SamplesToAdd = 1;
0263                   if (Dump >= 4)
0264                     cout << "            ......................................" << endl;
0265                   if (Dump >= 4)
0266                     for (int k = 0; k < theDigis->size(); k++) {
0267                       const HcalQIESample QIE = theDigis->sample(k);
0268                       if (k >= theDigis->presamples() && k < theDigis->presamples() + SamplesToAdd)
0269                         SumDigiCharge += QIE.nominal_fC();
0270                       if (k < theDigis->presamples() - 1)
0271                         EstimatedPedestal += QIE.nominal_fC() * SamplesToAdd / (theDigis->presamples() - 1);
0272                       cout << "            Digi: " << k << ", cap ID = " << QIE.capid()
0273                            << ": ADC Counts = " << QIE.adc() << ", nominal fC = " << QIE.nominal_fC() << endl;
0274                     }
0275                   if (Dump >= 4)
0276                     cout << "            1 Digi fC =" << SumDigiCharge << ", est. ped. fC=" << EstimatedPedestal
0277                          << ", est. GeV/fc=" << theRecHit->energy() / (SumDigiCharge - EstimatedPedestal) << endl;
0278                   if (Dump >= 4)
0279                     cout << "            ......................................" << endl;
0280                 }
0281               }
0282               if (ShowECal && DetNum == DetId::Ecal) {
0283                 int EcalNum = RecHitDetID.subdetId();
0284                 if (EcalNum == 1) {
0285                   EBDetId EcalID = RecHitDetID;
0286                   EBRecHitCollection::const_iterator theRecHit = EBRecHits->find(EcalID);
0287                   EBDigiCollection::const_iterator theDigis = EBDigis->find(EcalID);
0288                   sumRecHitE += theRecHit->energy();
0289                   cout << "         RecHit " << j << ": EB, ieta=" << EcalID.ieta() << ", iphi=" << EcalID.iphi()
0290                        << ", SM=" << EcalID.ism() << ", energy=" << theRecHit->energy()
0291                        << ", All Digis=" << theDigis->size() << endl;
0292                   if (Dump >= 4)
0293                     cout << "            ......................................" << endl;
0294                   if (Dump >= 4)
0295                     for (unsigned int k = 0; k < theDigis->size(); k++) {
0296                       EBDataFrame frame(*theDigis);
0297                       const EcalMGPASample MGPA = frame.sample(k);
0298                       cout << "            Digi: " << k << ": ADC Sample = " << MGPA.adc()
0299                            << ", Gain ID = " << MGPA.gainId() << endl;
0300                     }
0301                   if (Dump >= 4)
0302                     cout << "            ......................................" << endl;
0303                 } else if (EcalNum == 2) {
0304                   EEDetId EcalID = RecHitDetID;
0305                   EERecHitCollection::const_iterator theRecHit = EERecHits->find(EcalID);
0306                   EEDigiCollection::const_iterator theDigis = EEDigis->find(EcalID);
0307                   sumRecHitE += theRecHit->energy();
0308                   cout << "         RecHit " << j << ": EE, ix=" << EcalID.ix() << ", iy=" << EcalID.iy()
0309                        << ", energy=" << theRecHit->energy() << ", All Digis=" << theDigis->size() << endl;
0310                   if (Dump >= 4)
0311                     cout << "            ......................................" << endl;
0312                   if (Dump >= 4)
0313                     for (unsigned int k = 0; k < theDigis->size(); k++) {
0314                       EEDataFrame frame(*theDigis);
0315                       const EcalMGPASample MGPA = frame.sample(k);
0316                       cout << "            Digi: " << k << ": ADC Sample = " << MGPA.adc()
0317                            << ", Gain ID = " << MGPA.gainId() << endl;
0318                     }
0319                   if (Dump >= 4)
0320                     cout << "            ......................................" << endl;
0321                 }
0322               }
0323             }
0324           if (Dump >= 3) {
0325             if (abs(ietaTower) == 28 || abs(ietaTower) == 29) {
0326               cout << "         Splitted Sum of RecHit Energies=" << sumRecHitE
0327                    << ", CaloTower energy=" << theTower->energy() << endl;
0328             } else {
0329               cout << "         Sum of RecHit Energies=" << sumRecHitE << ", CaloTower energy=" << theTower->energy()
0330                    << endl;
0331             }
0332           }
0333           if (Dump >= 3)
0334             cout << "         ------------------------------------------------" << endl;
0335         }
0336       if (Dump >= 2)
0337         cout << "      Sum of tower energies=" << sumTowerE << ", CaloJet energy=" << jet->energy() << endl;
0338       jetInd++;
0339       if (Dump >= 2)
0340         cout << "      =====================================================" << endl;
0341     }
0342   evtCount++;
0343   if (Dump >= 1)
0344     cout << "   *********************************************************" << endl;
0345 }
0346 
0347 void JetToDigiDump::endJob() {}
0348 #include "FWCore/Framework/interface/MakerMacros.h"
0349 DEFINE_FWK_MODULE(JetToDigiDump);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!