Reference Manual

Go to the documentation of this file.
00001 // JetToDigiDump.cc
00002 // Description:  Prints out Jets, consituent CaloTowers, constituent RecHits and associated Digis (the digis for HCAL only).
00003 //               The user can specify which level in the config file:
00004 //               DumpLevel="Jets":    Printout of jets and their kinematic quantities.
00005 //               DumpLevel="Towers":  Nested Printout of jets and their constituent CaloTowers
00006 //               DumpLevel="RecHits": Nested Printout of jets, constituent CaloTowers and constituent RecHits
00007 //               DumpLevel="Digis":   Nested Printout of jets, constituent CaloTowers, RecHits and all the HCAL digis 
00008 //                                    associated with the RecHit channel (no links exist to go back to actual digis used).
00009 //               Does simple sanity checks on energy sums at each level: jets=sum of towers, tower=sum of RecHits.
00010 //               Does quick and dirty estimate of the fC/GeV factor that was applied to make the RecHit from the Digis.
00011 //               
00012 // Author: Robert M. Harris
00013 // Date:  19 - October - 2006
00014 // 
00015 #include "RecoJets/JetAnalyzers/interface/JetToDigiDump.h"
00016 #include "DataFormats/JetReco/interface/CaloJetCollection.h"
00017 #include "DataFormats/JetReco/interface/CaloJet.h"
00018 //in CaloJet: #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
00019 // just for the CaloTowerPtr declaration:
00020 // in CaloTowerCollection: #include "DataFormats/CaloTowers/interface/CaloTower.h"
00021 // in CaloTowerCollection: #include "DataFormats/CaloTowers/interface/CaloTowerFwd.h"
00022 #include "DataFormats/DetId/interface/DetId.h"
00023 #include "DataFormats/HcalDetId/interface/HcalDetId.h"
00024 #include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
00025 #include "DataFormats/HcalRecHit/interface/HBHERecHit.h"
00026 #include "DataFormats/EcalDetId/interface/EBDetId.h"
00027 #include "DataFormats/EcalDetId/interface/EEDetId.h"
00028 #include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
00029 #include "DataFormats/HcalDigi/interface/HcalDigiCollections.h"
00030 #include "DataFormats/EcalDigi/interface/EBDataFrame.h"
00031 #include "DataFormats/EcalDigi/interface/EcalDigiCollections.h"
00032 #include "FWCore/Framework/interface/Event.h"
00033 #include "FWCore/ParameterSet/interface/ParameterSet.h"
00034 #include <TROOT.h>
00035 #include <TSystem.h>
00036 #include <TFile.h>
00037 #include <TCanvas.h>
00038 #include <cmath>
00039 using namespace edm;
00040 using namespace reco;
00041 using namespace std;
00042 
00043 JetToDigiDump::JetToDigiDump( const ParameterSet & cfg ) :
00044   DumpLevel( cfg.getParameter<string>( "DumpLevel" ) ),
00045   CaloJetAlg( cfg.getParameter<string>( "CaloJetAlg" ) ),
00046   DebugLevel( cfg.getParameter<int>( "DebugLevel" ) ),
00047   ShowECal( cfg.getParameter<bool>( "ShowECal" ) )
00048   {
00049 }
00050 
00051 void JetToDigiDump::beginJob( ) {
00052   if(DumpLevel=="Jets")
00053   {
00054     cout << "Dump of Jets" << endl;
00055     Dump=1;
00056   }
00057   else if(DumpLevel=="Towers")
00058   {
00059     cout << "Dump of Jets and constituent CaloTowers" << endl;    
00060     Dump=2;
00061   }
00062   else if(DumpLevel=="RecHits")
00063   {
00064     cout << "Dump of Jets, constituent CaloTowers, and constituent RecHits" << endl;    
00065     Dump=3;
00066   }
00067   else if(DumpLevel=="Digis")
00068   {
00069     cout << "Dump of Jets, constituent CaloTowers, constituent RecHits and associated Digis" << endl;    
00070     Dump=4;
00071   }
00072   cout << "Jet Algorithm being dumped is " << CaloJetAlg << endl;
00073   cout<<"Debug level is " << DebugLevel << endl;
00074   //Initialize some stuff
00075   evtCount = 0;
00076 }
00077 
00078 void JetToDigiDump::analyze( const Event& evt, const EventSetup& es ) {
00079 
00080   int jetInd;
00081   Handle<CaloJetCollection> caloJets;
00082   Handle<CaloTowerCollection> caloTowers;
00083   Handle<HBHERecHitCollection> HBHERecHits;
00084   Handle<HORecHitCollection> HORecHits;
00085   Handle<HFRecHitCollection> HFRecHits;
00086   Handle<EBRecHitCollection> EBRecHits;
00087   Handle<EERecHitCollection> EERecHits;
00088   Handle<HBHEDigiCollection> HBHEDigis;
00089   Handle<HODigiCollection> HODigis;
00090   Handle<HFDigiCollection> HFDigis;
00091   Handle<EEDigiCollection> EEDigis;
00092   Handle<EBDigiCollection> EBDigis;
00093   // Old:
00094   //Handle<edm::SortedCollection<EBDataFrame> > EBDigis;
00095   //Handle<edm::SortedCollection<EBDataFrame> > EEDigis;
00096    
00097   //Find the CaloTowers in leading CaloJets
00098   if (DebugLevel) cout<<"Getting caloJets"<<endl;
00099 
00100   evt.getByLabel( CaloJetAlg, caloJets );
00101   if (Dump >= 2) evt.getByLabel( "towerMaker", caloTowers );
00102   if (Dump >= 3) {
00103     if (DebugLevel) cout<<"Getting recHits"<<endl;
00104     evt.getByLabel( "hbhereco", HBHERecHits );
00105     evt.getByLabel( "horeco", HORecHits );
00106     evt.getByLabel( "hfreco", HFRecHits );
00107     evt.getByLabel( "ecalRecHit", "EcalRecHitsEB", EBRecHits );
00108     evt.getByLabel( "ecalRecHit", "EcalRecHitsEE", EERecHits );
00109     if (DebugLevel) cout<<"# of hits gotten - HBHE: "<<HBHERecHits->size()<<endl;
00110     evt.getByLabel( "hcalDigis", HBHEDigis );
00111     evt.getByLabel( "hcalDigis", HODigis );
00112     evt.getByLabel( "hcalDigis", HFDigis );
00113     if (DebugLevel) cout<<"# of digis gotten - HBHE: "<<HBHEDigis->size()<<endl;
00114     if (ShowECal) {
00115       evt.getByLabel( "ecalDigis", "ebDigis", EBDigis );
00116       evt.getByLabel( "ecalDigis", "eeDigis", EEDigis );
00117     }
00118   }
00119     
00120   cout << endl << "Evt: "<<evtCount <<", Num Jets=" <<caloJets->end() - caloJets->begin() << endl;
00121   if(Dump>=1)cout <<"   *********************************************************" <<endl;
00122   jetInd = 0;
00123   if(Dump>=1)for( CaloJetCollection::const_iterator jet = caloJets->begin(); jet != caloJets->end(); ++ jet ) {
00124     //2_1_?    std::vector <CaloTowerPtr> towers = jet->getCaloConstituents ();
00125     //2_0_7"
00126     std::vector <CaloTowerPtr> towers = jet->getCaloConstituents ();
00127     int nConstituents= towers.size();
00128     cout <<"   Jet: "<<jetInd<<", eta="<<jet->eta()<<", phi="<<jet->phi()<<", pt="<<jet->pt()<<\
00129     ",E="<<jet->energy()<<", EB E="<<jet->emEnergyInEB()<<" ,HB E="<<jet->hadEnergyInHB()<<\
00130     ", HO E="<<jet->hadEnergyInHO()<<" ,EE E="<< jet->emEnergyInEE()\
00131      <<", HE E="<<jet->hadEnergyInHE()<<", HF E="<<jet->hadEnergyInHF()+jet->emEnergyInHF()<<", Num Towers="<<nConstituents<<endl;
00132     if(Dump>=2)cout <<"      ====================================================="<<endl;
00133     float sumTowerE = 0.0;
00134     if(Dump>=2)for (int i = 0; i <nConstituents ; i++) {
00135        CaloTowerCollection::const_iterator theTower=caloTowers->find(towers[i]->id());  //Find the tower from its CaloTowerDetID        
00136        if (theTower == caloTowers->end()) {cerr<<"Bug? Can't find the tower"<<endl; return;}
00137        int ietaTower = towers[i]->id().ieta();
00138        int iphiTower = towers[i]->id().iphi();
00139        sumTowerE += theTower->energy();
00140        size_t numRecHits = theTower->constituentsSize();
00141        cout << "      Tower " << i <<": ieta=" << ietaTower <<  ", eta=" << theTower->eta() <<", iphi=" << iphiTower << ", phi=" << theTower->phi() << \
00142        ", energy=" << theTower->energy() << ", EM=" << theTower->emEnergy()<< ", HAD=" << theTower->hadEnergy()\
00143        << ", HO=" << theTower->outerEnergy() <<",  Num Rec Hits =" << numRecHits << endl;
00144        if(Dump>=3)cout << "         ------------------------------------------------"<<endl;
00145        float sumRecHitE = 0.0;
00146        if(Dump>=3)for(size_t j = 0; j <numRecHits ; j++) {
00147           DetId RecHitDetID=theTower->constituent(j);
00148           DetId::Detector DetNum=RecHitDetID.det();
00149           if( DetNum == DetId::Hcal ){
00150             //cout << "RecHit " << j << ": Detector = " << DetNum << ": Hcal " << endl;
00151             HcalDetId HcalID = RecHitDetID;
00152             HcalSubdetector HcalNum = HcalID.subdet();
00153             if(  HcalNum == HcalBarrel ){
00154               HBHERecHitCollection::const_iterator theRecHit=HBHERecHits->find(HcalID);     
00155               sumRecHitE += theRecHit->energy();
00156               HBHEDigiCollection::const_iterator theDigis=HBHEDigis->find(HcalID);
00157               cout << "         RecHit: " << j << ": HB, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()<<      
00158               ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy() << ", time=" <<\
00159               theRecHit->time() <<", All Digis=" << theDigis->size() << ", presamples =" <<\
00160               theDigis->presamples() <<endl;              
00161               /* const HcalElectronicsId HW_ID = theDigis->elecId();
00162                 cout << "Digi: Index=" << HW_ID.linearIndex() << ", raw ID=" <<  HW_ID.rawId() << ", fiberChan=" << HW_ID.fiberChanId() <<  ", fiberInd=" <<  HW_ID.fiberIndex() \
00163                      << ", HTR chan=" <<  HW_ID.htrChanId() << ", HTR Slot=" <<   HW_ID.htrSlot() << ", HDR top/bot=" << HW_ID.htrTopBottom() \
00164                 << ", VME crate=" <<  HW_ID.readoutVMECrateId() << ", DCC=" << HW_ID.dccid() << ", DCC spigot=" <<  HW_ID.spigot() << endl;
00165               */
00166               float SumDigiCharge = 0.0;
00167               float EstimatedPedestal=0.0;
00168               int SamplesToAdd = 4;
00169               if(Dump>=4)cout << "            ......................................"<<endl;
00170               if(Dump>=4)for(int k=0; k<theDigis->size(); k++){
00171                 const HcalQIESample QIE = theDigis->sample(k);
00172                 if(k>=theDigis->presamples()&&k<theDigis->presamples()+SamplesToAdd)SumDigiCharge+=QIE.nominal_fC();
00173                 if(k<theDigis->presamples()-1)EstimatedPedestal+=QIE.nominal_fC()*SamplesToAdd/(theDigis->presamples()-1);
00174                 cout << "            Digi: " << k <<  ", cap ID = " << QIE.capid() << ": ADC Counts = " << QIE.adc() <<  ", nominal fC = " << QIE.nominal_fC() <<endl;
00175               }
00176               if(Dump>=4)cout << "            4 Digi fC ="<<SumDigiCharge<<", est. ped. fC="<<EstimatedPedestal<<", est. GeV/fc="<<theRecHit->energy()/(SumDigiCharge-EstimatedPedestal) << endl;
00177               if(Dump>=4)cout << "            ......................................"<<endl;
00178             }
00179             else if(  HcalNum == HcalEndcap  ){
00180               HBHERecHitCollection::const_iterator theRecHit=HBHERecHits->find(HcalID);     
00181               if( (abs(HcalID.ieta())==28||abs(HcalID.ieta())==29)&&HcalID.depth()==3){
00182                 sumRecHitE += theRecHit->energy()/2;  //Depth 3 split over tower 28 & 29
00183               }
00184               else{
00185                 sumRecHitE += theRecHit->energy();
00186               }
00187               HBHEDigiCollection::const_iterator theDigis=HBHEDigis->find(HcalID);
00188               cout << "         RecHit: " << j << ": HE, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()<<      
00189               ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy() << ", time=" <<\
00190               theRecHit->time() <<", All Digis=" << theDigis->size() << ", presamples =" <<\
00191               theDigis->presamples() <<endl;              
00192               float SumDigiCharge = 0.0;
00193               float EstimatedPedestal=0.0;
00194               int SamplesToAdd = 4;
00195               if(Dump>=4)cout << "            ......................................"<<endl;
00196               if(Dump>=4)for(int k=0; k<theDigis->size(); k++){
00197                 const HcalQIESample QIE = theDigis->sample(k);
00198                 if(k>=theDigis->presamples()&&k<theDigis->presamples()+SamplesToAdd)SumDigiCharge+=QIE.nominal_fC();
00199                 if(k<theDigis->presamples()-1)EstimatedPedestal+=QIE.nominal_fC()*SamplesToAdd/(theDigis->presamples()-1);
00200                 cout << "            Digi: " << k <<  ", cap ID = " << QIE.capid() << ": ADC Counts = " << QIE.adc() <<  ", nominal fC = " << QIE.nominal_fC() <<endl;
00201               }
00202               if(Dump>=4)cout << "            4 Digi fC ="<<SumDigiCharge<<", est. ped. fC="<<EstimatedPedestal<<", est. GeV/fc="<<theRecHit->energy()/(SumDigiCharge-EstimatedPedestal) << endl;
00203               if(Dump>=4)cout << "            ......................................"<<endl;
00204             }        
00205             else if(  HcalNum == HcalOuter  ){
00206               HORecHitCollection::const_iterator theRecHit=HORecHits->find(HcalID);         
00207               sumRecHitE += theRecHit->energy();
00208               HODigiCollection::const_iterator theDigis=HODigis->find(HcalID);
00209               cout << "         RecHit: " << j << ": HO, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()<<      
00210               ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy() << ", time=" <<\
00211               theRecHit->time() <<", All Digis=" << theDigis->size() << ", presamples =" <<\
00212               theDigis->presamples() <<endl;              
00213               float SumDigiCharge = 0.0;
00214               float EstimatedPedestal=0.0;
00215               int SamplesToAdd = 4;
00216               if(Dump>=4)cout << "            ......................................"<<endl;
00217               if(Dump>=4)for(int k=0; k<theDigis->size(); k++){
00218                 const HcalQIESample QIE = theDigis->sample(k);
00219                 if(k>=theDigis->presamples()&&k<theDigis->presamples()+SamplesToAdd)SumDigiCharge+=QIE.nominal_fC();
00220                 if(k<theDigis->presamples()-1)EstimatedPedestal+=QIE.nominal_fC()*SamplesToAdd/(theDigis->presamples()-1);
00221                 cout << "            Digi: " << k <<  ", cap ID = " << QIE.capid() << ": ADC Counts = " << QIE.adc() <<  ", nominal fC = " << QIE.nominal_fC() <<endl;
00222               }
00223               if(Dump>=4)cout << "            4 Digi fC ="<<SumDigiCharge<<", est. ped. fC="<<EstimatedPedestal<<", est. GeV/fc="<<theRecHit->energy()/(SumDigiCharge-EstimatedPedestal) << endl;
00224               if(Dump>=4)cout << "            ......................................"<<endl;
00225             }        
00226             else if(  HcalNum == HcalForward  ){
00227               HFRecHitCollection::const_iterator theRecHit=HFRecHits->find(HcalID);         
00228               sumRecHitE += theRecHit->energy();
00229               HFDigiCollection::const_iterator theDigis=HFDigis->find(HcalID);
00230               cout << "         RecHit: " << j << ": HF, ieta=" << HcalID.ieta() << ", iphi=" << HcalID.iphi()<<      
00231               ", depth=" << HcalID.depth() << ", energy=" << theRecHit->energy() << ", time=" <<\
00232               theRecHit->time() <<", All Digis=" << theDigis->size() << ", presamples =" <<\
00233               theDigis->presamples() <<endl;              
00234               float SumDigiCharge = 0.0;
00235               float EstimatedPedestal=0.0;            
00236               int SamplesToAdd = 1;
00237               if(Dump>=4)cout << "            ......................................"<<endl;
00238               if(Dump>=4)for(int k=0; k<theDigis->size(); k++){
00239                 const HcalQIESample QIE = theDigis->sample(k);
00240                 if(k>=theDigis->presamples()&&k<theDigis->presamples()+SamplesToAdd)SumDigiCharge+=QIE.nominal_fC();
00241                 if(k<theDigis->presamples()-1)EstimatedPedestal+=QIE.nominal_fC()*SamplesToAdd/(theDigis->presamples()-1);
00242                 cout << "            Digi: " << k <<  ", cap ID = " << QIE.capid() << ": ADC Counts = " << QIE.adc() <<  ", nominal fC = " << QIE.nominal_fC() <<endl;
00243               }
00244               if(Dump>=4)cout << "            1 Digi fC ="<<SumDigiCharge<<", est. ped. fC="<<EstimatedPedestal<<", est. GeV/fc="<<theRecHit->energy()/(SumDigiCharge-EstimatedPedestal) << endl;
00245               if(Dump>=4)cout << "            ......................................"<<endl;
00246            }                                  
00247           }
00248           if( ShowECal && DetNum == DetId::Ecal ){
00249             int EcalNum =  RecHitDetID.subdetId();
00250             if( EcalNum == 1 ){
00251               EBDetId EcalID = RecHitDetID;
00252               EBRecHitCollection::const_iterator theRecHit=EBRecHits->find(EcalID);         
00253               EBDigiCollection::const_iterator theDigis=EBDigis->find(EcalID);
00254               sumRecHitE += theRecHit->energy();
00255               cout << "         RecHit " << j << ": EB, ieta=" << EcalID.ieta() <<  ", iphi=" << EcalID.iphi() <<  ", SM=" << EcalID.ism() << ", energy=" << theRecHit->energy() <<", All Digis=" << theDigis->size()<< endl;
00256               if(Dump>=4)cout << "            ......................................"<<endl;
00257               if(Dump>=4)for(unsigned int k=0; k<theDigis->size(); k++){
00258                 EBDataFrame frame (*theDigis);
00259                 const EcalMGPASample MGPA = frame.sample(k);
00260                 cout << "            Digi: " << k <<   ": ADC Sample = " << MGPA.adc() << ", Gain ID = "<< MGPA.gainId() <<endl;
00261               }
00262               if(Dump>=4)cout << "            ......................................"<<endl;
00263             }
00264             else if(  EcalNum == 2 ){
00265               EEDetId EcalID = RecHitDetID;
00266               EERecHitCollection::const_iterator theRecHit=EERecHits->find(EcalID);         
00267               EEDigiCollection::const_iterator theDigis=EEDigis->find(EcalID);
00268               sumRecHitE += theRecHit->energy();
00269               cout << "         RecHit " << j << ": EE, ix=" << EcalID.ix() <<  ", iy=" << EcalID.iy() << ", energy=" << theRecHit->energy() << ", All Digis=" << theDigis->size()<< endl;
00270               if(Dump>=4)cout << "            ......................................"<<endl;
00271               if(Dump>=4)for(unsigned int k=0; k<theDigis->size(); k++){
00272                 EEDataFrame frame (*theDigis);
00273                 const EcalMGPASample MGPA = frame.sample(k);
00274                 cout << "            Digi: " << k <<   ": ADC Sample = " << MGPA.adc() << ", Gain ID = "<< MGPA.gainId() <<endl;
00275               }
00276               if(Dump>=4)cout << "            ......................................"<<endl;
00277             }
00278           }
00279        }
00280        if(Dump>=3){
00281          if( abs(ietaTower)==28||abs(ietaTower)==29){
00282              cout << "         Splitted Sum of RecHit Energies=" << sumRecHitE <<", CaloTower energy=" << theTower->energy() <<  endl;
00283          }
00284          else{
00285              cout << "         Sum of RecHit Energies=" << sumRecHitE <<", CaloTower energy=" << theTower->energy() <<  endl;
00286          }
00287        }
00288        if(Dump>=3)cout << "         ------------------------------------------------"<<endl;
00289     }
00290     if(Dump>=2)cout << "      Sum of tower energies=" << sumTowerE << ", CaloJet energy=" << jet->energy() <<  endl;
00291     jetInd++;
00292     if(Dump>=2)cout <<"      ====================================================="<<endl;
00293   }
00294   evtCount++;    
00295   if(Dump>=1)cout <<"   *********************************************************" <<endl;
00296 
00297 }
00298 
00299 void JetToDigiDump::endJob() {
00300 
00301 
00302 }
00303 #include "FWCore/Framework/interface/MakerMacros.h"
00304 DEFINE_FWK_MODULE(JetToDigiDump);
/data/refman/pasoursint/CMSSW_6_1_2_SLHC2_patch1/src/RecoJets/JetAnalyzers/src/JetToDigiDump.cc
