JetAnalysis/src/RecHitComparison.cc

// -*- C++ -*-
//
// Package:    RecHitComparison
// Class:      RecHitComparison
// 
/**\class RecHitComparison RecHitComparison.cc CmsHi/RecHitComparison/src/RecHitComparison.cc

 Description: [one line class summary]

 Implementation:
     [Notes on implementation]
*/
//
// Original Author:  Yetkin Yilmaz
//         Created:  Tue Sep  7 11:38:19 EDT 2010
// $Id: RecHitComparison.cc,v 1.9 2011/03/30 12:13:22 yilmaz Exp $
//
//


// system include files
#include <memory>
#include <vector>
#include <iostream>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/DetId/interface/DetId.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"

#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "DataFormats/Math/interface/deltaR.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/FWLite/interface/ChainEvent.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
#include "DataFormats/CaloTowers/interface/CaloTower.h"
#include "DataFormats/HeavyIonEvent/interface/CentralityProvider.h"
#include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"

#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"

#include "TNtuple.h"

using namespace std;

//
// class declaration
//

class RecHitComparison : public edm::EDAnalyzer {
   public:
      explicit RecHitComparison(const edm::ParameterSet&);
      ~RecHitComparison();


   private:
      virtual void beginJob() ;
      virtual void analyze(const edm::Event&, const edm::EventSetup&);
      virtual void endJob() ;

      // ----------member data ---------------------------
   edm::Handle<vector<double> > ktRhos;
   edm::Handle<vector<double> > akRhos;

  edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits1;
  edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits2;
  edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits1;
  edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits2;

  edm::Handle<HFRecHitCollection> hfHits1;
  edm::Handle<HFRecHitCollection> hfHits2;
  edm::Handle<HBHERecHitCollection> hbheHits1;
  edm::Handle<HBHERecHitCollection> hbheHits2;

   edm::Handle<reco::BasicClusterCollection> bClusters1;
   edm::Handle<reco::BasicClusterCollection> bClusters2;


   typedef vector<EcalRecHit>::const_iterator EcalIterator;
  typedef vector<HFRecHit>::const_iterator HFIterator;
  typedef vector<HBHERecHit>::const_iterator HBHEIterator;

   edm::Handle<reco::CaloJetCollection> signalJets;

  edm::InputTag HcalRecHitHFSrc1_;
  edm::InputTag HcalRecHitHFSrc2_;
  edm::InputTag HcalRecHitHBHESrc1_;
  edm::InputTag HcalRecHitHBHESrc2_;
  edm::InputTag EBSrc1_;
  edm::InputTag EBSrc2_;
  edm::InputTag EESrc1_;
  edm::InputTag EESrc2_;

   edm::InputTag BCSrc1_;
   edm::InputTag BCSrc2_;

   edm::InputTag signalTag_;

   TNtuple* ntBC;
   TNtuple* ntEB;
  TNtuple* ntEE;
  TNtuple* ntHBHE;
  TNtuple* ntHF;
   TNtuple* ntjet;

   double cone;
   bool jetsOnly_;
   bool doBasicClusters_;
  bool doJetCone_;

   edm::Service<TFileService> fs;
   CentralityProvider* centrality_;
   const CaloGeometry *geo;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
RecHitComparison::RecHitComparison(const edm::ParameterSet& iConfig) :
   centrality_(0),
   geo(0),
   cone(0.5)
{
   //now do what ever initialization is needed
   jetsOnly_ = iConfig.getUntrackedParameter<bool>("jetsOnly",false);
   doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);

   doJetCone_ = iConfig.getUntrackedParameter<bool>("doJetCone",false);
   signalTag_ = iConfig.getUntrackedParameter<edm::InputTag>("signalJets",edm::InputTag("iterativeCone5CaloJets","","SIGNAL"));

   if(!doJetCone_) jetsOnly_ = 0;

  HcalRecHitHFSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc1",edm::InputTag("hfreco"));
  HcalRecHitHFSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc2",edm::InputTag("hfreco"));
  HcalRecHitHBHESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc1",edm::InputTag("hbhereco"));
  HcalRecHitHBHESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc2",edm::InputTag("hbhereco"));
  EBSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECOBKG"));
  EBSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","S"));
  EESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEE","RECO"));
  EESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEE","SIGNALONLY"));
  BCSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
  BCSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));

}


RecHitComparison::~RecHitComparison()
{
 
   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
RecHitComparison::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
{
   if(!centrality_) centrality_ = new CentralityProvider(iSetup);
   if(!geo){
      edm::ESHandle<CaloGeometry> pGeo;
      iSetup.get<CaloGeometryRecord>().get(pGeo);
      geo = pGeo.product();
   }

   using namespace edm;
   ev.getByLabel(EBSrc1_,ebHits1);
   ev.getByLabel(EBSrc2_,ebHits2);

   if(doJetCone_) ev.getByLabel(signalTag_,signalJets);

   ev.getByLabel(HcalRecHitHFSrc1_,hfHits1);
   ev.getByLabel(HcalRecHitHFSrc2_,hfHits2);
   ev.getByLabel(HcalRecHitHBHESrc1_,hbheHits1);
   ev.getByLabel(HcalRecHitHBHESrc2_,hbheHits2);
   ev.getByLabel(EESrc1_,eeHits1);
   ev.getByLabel(EESrc2_,eeHits2);

   if(doBasicClusters_){
      ev.getByLabel(BCSrc1_,bClusters1);
      ev.getByLabel(BCSrc2_,bClusters2);
   }

   centrality_->newEvent(ev,iSetup);

   double hf = centrality_->centralityValue();
   int bin = centrality_->getBin();

   vector<double> fFull;
   vector<double> f05;
   vector<double> f1;
   vector<double> f15;
   vector<double> f2;
   vector<double> f25;
   vector<double> f3;

   int njets = 0;

   if(doJetCone_) njets = signalJets->size();
   fFull.reserve(njets);
   f05.reserve(njets);
   f1.reserve(njets);
   f15.reserve(njets);
   f2.reserve(njets);
   f25.reserve(njets);
   f3.reserve(njets);

   if(doJetCone_){
     for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
       fFull.push_back(0);
       f05.push_back(0);
       f1.push_back(0);
       f15.push_back(0);
       f2.push_back(0);
       f25.push_back(0);
       f3.push_back(0);
     }
   }

   for(unsigned int j1 = 0 ; j1 < ebHits1->size(); ++j1){

      const EcalRecHit& jet1 = (*ebHits1)[j1];
      double e1 = jet1.energy();

      const GlobalPoint& pos1=geo->getPosition(jet1.id());
      double eta1 = pos1.eta();
      double phi1 = pos1.phi();
      double et1 = e1*sin(pos1.theta());

      double drjet = -1;
      double jetpt = -1;
      bool isjet = false;
      int matchedJet = -1;

      if(doJetCone_){
        for(unsigned int j = 0 ; j < signalJets->size(); ++j){
          const reco::CaloJet & jet = (*signalJets)[j];
          double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
          if(dr < cone){
            jetpt = jet.pt();
            drjet = dr;
            isjet = true;
            matchedJet = j;
            fFull[j] += et1;
            
            if(et1 > 0.5){
              f05[j] += et1;
            }
            if(et1 > 1.){
               f1[j] += et1;
            }
            if(et1 > 1.5){
               f15[j] += et1;
            }
            if(et1 > 2){
               f2[j] += et1;
            }
            if(et1 > 2.5){
               f25[j] += et1;
            }
            if(et1 > 3.){
               f3[j] += et1;
            }
          }
        }
      }   
      
      GlobalPoint pos2;
      double e2 = -1;
      EcalIterator hitit = ebHits2->find(jet1.id());
      if(hitit != ebHits2->end()){
         e2 = hitit->energy();
         pos2=geo->getPosition(hitit->id());
      }else{
         e2 = 0;
         pos2 = pos1; 
      }

      double eta2 = pos2.eta();
      double phi2 = pos2.eta();
      double et2 = e2*sin(pos2.theta());
      if(!jetsOnly_ ||  isjet) ntEB->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
   }

   for(unsigned int i = 0; i < eeHits1->size(); ++i){
     const EcalRecHit & jet1= (*eeHits1)[i];
     double e1 = jet1.energy();
     const GlobalPoint& pos1=geo->getPosition(jet1.id());
     double eta1 = pos1.eta();
     double phi1 = pos1.phi();
     double et1 = e1*sin(pos1.theta());
     double drjet = -1;
     double jetpt = -1;
     bool isjet = false;
     int matchedJet = -1;
     if(doJetCone_){
       for(unsigned int j = 0 ; j < signalJets->size(); ++j){
         const reco::CaloJet & jet = (*signalJets)[j];
         double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
         if(dr < cone){
           jetpt = jet.pt();
           drjet = dr;
           isjet = true;
           matchedJet = j;
         }
       }
     }
     
     GlobalPoint pos2;
     double e2 = -1;
     EcalIterator hitit = eeHits2->find(jet1.id());
     if(hitit != eeHits2->end()){
       e2 = hitit->energy();
       pos2=geo->getPosition(hitit->id());
     }else{
       e2 = 0;
       pos2 = pos1;
     }
     double eta2 = pos2.eta();
     double phi2 = pos2.eta();
     double et2 = e2*sin(pos2.theta());
     if(!jetsOnly_ || isjet) ntEE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
   }
   
   for(unsigned int i = 0; i < hbheHits1->size(); ++i){
     const HBHERecHit & jet1= (*hbheHits1)[i];
     double e1 = jet1.energy();
     const GlobalPoint& pos1=geo->getPosition(jet1.id());
     double eta1 = pos1.eta();
     double phi1 = pos1.phi();
     double et1 = e1*sin(pos1.theta());
     double drjet = -1;
     double jetpt = -1;
     bool isjet = false;
     int matchedJet = -1;
     if(doJetCone_){
       for(unsigned int j = 0 ; j < signalJets->size(); ++j){
         const reco::CaloJet & jet = (*signalJets)[j];
         double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
         if(dr < cone){
           jetpt = jet.pt();
           drjet = dr;
           isjet = true;
           matchedJet = j;
         }
       }
     }

     GlobalPoint pos2;
     double e2 = -1;
     HBHEIterator hitit = hbheHits2->find(jet1.id());
     if(hitit != hbheHits2->end()){
       e2 = hitit->energy();
       pos2=geo->getPosition(hitit->id());
     }else{
       e2 = 0;
       pos2 = pos1;
     }
     double eta2 = pos2.eta();
     double phi2 = pos2.eta();
     double et2 = e2*sin(pos2.theta());
     if(!jetsOnly_ || isjet) ntHBHE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
   }

   for(unsigned int i = 0; i < hfHits1->size(); ++i){
     const HFRecHit & jet1= (*hfHits1)[i];
     double e1 = jet1.energy();
     const GlobalPoint& pos1=geo->getPosition(jet1.id());
     double eta1 = pos1.eta();
     double phi1 = pos1.phi();
     double et1 = e1*sin(pos1.theta());
     double drjet = -1;
     double jetpt = -1;
     bool isjet = false;
     int matchedJet = -1;
     if(doJetCone_){
       for(unsigned int j = 0 ; j < signalJets->size(); ++j){
         const reco::CaloJet & jet = (*signalJets)[j];
         double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
         if(dr < cone){
           jetpt = jet.pt();
           drjet = dr;
           isjet = true;
           matchedJet = j;
         }
       }
     }
     GlobalPoint pos2;
     double e2 = -1;
     HFIterator hitit = hfHits2->find(jet1.id());
     if(hitit != hfHits2->end()){
       e2 = hitit->energy();
       pos2=geo->getPosition(hitit->id());
     }else{
       e2 = 0;
       pos2 = pos1;
     }
     double eta2 = pos2.eta();
     double phi2 = pos2.eta();
     double et2 = e2*sin(pos2.theta());
     if(!jetsOnly_ || isjet) ntHF->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
   }

   if(doJetCone_){
     for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
       const reco::CaloJet & jet = (*signalJets)[j1];
       double em = (jet.emEnergyInEB() + jet.emEnergyInEE()) * sin(jet.theta());
       double emf = jet.emEnergyFraction(); 
       double pt = jet.pt();
       double eta = jet.eta();
       ntjet->Fill(bin,pt,eta,fFull[j1],f05[j1],f1[j1],f15[j1],f2[j1],f25[j1],f3[j1],em,emf);
     }
   }

}


// ------------ method called once each job just before starting event loop  ------------
void 
RecHitComparison::beginJob()
{
   ntEB = fs->make<TNtuple>("ntEB","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
   ntEE = fs->make<TNtuple>("ntEE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
   ntHBHE = fs->make<TNtuple>("ntHBHE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
   ntHF = fs->make<TNtuple>("ntHF","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");

   ntBC = fs->make<TNtuple>("ntBC","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");

   ntjet = fs->make<TNtuple>("ntjet","","bin:pt:eta:ethit:f05:f1:f15:f2:f25:f3:em:emf");
   
}

// ------------ method called once each job just after ending the event loop  ------------
void 
RecHitComparison::endJob() {
}

//define this as a plug-in
DEFINE_FWK_MODULE(RecHitComparison);
Revision:	1.10
Committed:	Wed Mar 30 21:21:10 2011 UTC (14 years, 1 month ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HiForest_V02_85, HiForest_V02_84, HiForest_V02_83, HiForest_V02_82, HiForest_V02_81, HiForest_V02_80, HiForest_V02_79, HiForest_V02_78, HiForest_V02_77, HiForest_V02_76, HiForest_V02_75, HiForest_V02_74, HiForest_V02_73, HiForest_V02_72, HiForest_V02_71, HiForest_V02_70, HiForest_V02_69, HiForest_V02_68, HiForest_V02_67, HiForest_V02_66, HiForest_V02_65, HiForest_V02_64, HiForest_V02_63, HiForest_V02_62, HiForest_V02_61, HiForest_V02_60, HiForest_V02_59, HiForest_V02_58, HiForest_V02_57, HiForest_V02_56, HiForest_V02_55, HiForest_V02_54, HiForest_V02_53, HiForest_V02_52, HiForest_V02_51, HiForest_V02_50, HiForest_V02_49, HiForest_V02_48, HiForest_V02_47, HiForest_V02_46, HiForest_V02_45, HiForest_V02_44, HiForest_V02_43, HiForest_V02_42, HiForest_V02_41, HiForest_V02_40, HiForest_V02_39, HiForest_V02_38, HiForest_V02_37, HiForest_V02_36, HiForest_V02_35, HiForest_V02_34, HiForest_V02_33, HiForest_V02_32, HiForest_V02_31, HiForest_V02_30, HiForest_V02_27, HiForest_V02_26, QM_2012, HiForest_V02_25, HiForest_V02_24, HiForest_V02_23, HiForest_V02_22, HiForest_V02_21, HiForest_V02_20, HiForest_V02_19, HiForest_V02_18, HiForest_V02_17, HiForest_V02_16, HiForest_V02_15, HiForest_V02_14, HiForest_V02_13, HiForest_V02_12, HiForest_V02_11, HiForest_V02_10, HiForest_V02_09, HiForest_V02_08, HiForest_V02_07, HiForest_V02_06, HiForest_V02_05, HiForest_V02_04, HiForest_V02_03, HiForest_V02_02, HiForest_V02_01, HiForest_V02_00, hi44X_02, hi413_03, hi441_1, hi441_0, hi413_11, hi413_10, hi413_09, hi413_08, hi413_07, hi413_06, hi413_05, hi413_04, hi413_02, hi39X_01, tag_d20110915, cmssw39x_base, cmssw39X_base, HEAD
Branch point for:	branch_44x, cmssw39x_branch
Changes since 1.9:	+26 -20 lines
Log Message:	jet optional
#	User	Rev	Content
1	yilmaz	1.1	// -- C++ --
2			//
3			// Package: RecHitComparison
4			// Class: RecHitComparison
5			//
6			/**\class RecHitComparison RecHitComparison.cc CmsHi/RecHitComparison/src/RecHitComparison.cc
7
8			Description: [one line class summary]
9
10			Implementation:
11			[Notes on implementation]
12			*/
13			//
14			// Original Author: Yetkin Yilmaz
15			// Created: Tue Sep 7 11:38:19 EDT 2010
16	yilmaz	1.10	// $Id: RecHitComparison.cc,v 1.9 2011/03/30 12:13:22 yilmaz Exp $
17	yilmaz	1.1	//
18			//
19
20
21			// system include files
22			#include <memory>
23			#include <vector>
24			#include <iostream>
25
26			// user include files
27			#include "FWCore/Framework/interface/Frameworkfwd.h"
28			#include "FWCore/Framework/interface/EDAnalyzer.h"
29
30			#include "FWCore/Framework/interface/Event.h"
31			#include "FWCore/Framework/interface/MakerMacros.h"
32			#include "FWCore/Framework/interface/ESHandle.h"
33			#include "FWCore/ParameterSet/interface/ParameterSet.h"
34
35			#include "DataFormats/DetId/interface/DetId.h"
36			#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
37			#include "Geometry/Records/interface/IdealGeometryRecord.h"
38			#include "Geometry/Records/interface/CaloGeometryRecord.h"
39			#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
40
41			#include "CommonTools/UtilAlgos/interface/TFileService.h"
42			#include "FWCore/ServiceRegistry/interface/Service.h"
43
44			#include "DataFormats/Math/interface/deltaR.h"
45
46			#include "DataFormats/Common/interface/Handle.h"
47			#include "DataFormats/FWLite/interface/ChainEvent.h"
48			#include "FWCore/Utilities/interface/InputTag.h"
49			#include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
50			#include "DataFormats/CaloTowers/interface/CaloTower.h"
51	yilmaz	1.6	#include "DataFormats/HeavyIonEvent/interface/CentralityProvider.h"
52	yilmaz	1.1	#include "SimDataFormats/HiGenData/interface/GenHIEvent.h"
53			#include "DataFormats/PatCandidates/interface/Jet.h"
54			#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
55	yilmaz	1.2	#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
56			#include "DataFormats/HcalDetId/interface/HcalDetId.h"
57	yilmaz	1.1
58	yilmaz	1.7	#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
59
60	yilmaz	1.1	#include "TNtuple.h"
61
62			using namespace std;
63
64			//
65			// class declaration
66			//
67
68			class RecHitComparison : public edm::EDAnalyzer {
69			public:
70			explicit RecHitComparison(const edm::ParameterSet&);
71			~RecHitComparison();
72
73
74			private:
75			virtual void beginJob() ;
76			virtual void analyze(const edm::Event&, const edm::EventSetup&);
77			virtual void endJob() ;
78
79			// ----------member data ---------------------------
80			edm::Handle<vector<double> > ktRhos;
81			edm::Handle<vector<double> > akRhos;
82
83	yilmaz	1.2	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits1;
84			edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits2;
85			edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits1;
86			edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits2;
87
88			edm::Handle<HFRecHitCollection> hfHits1;
89			edm::Handle<HFRecHitCollection> hfHits2;
90			edm::Handle<HBHERecHitCollection> hbheHits1;
91			edm::Handle<HBHERecHitCollection> hbheHits2;
92	yilmaz	1.1
93	yilmaz	1.7	edm::Handle<reco::BasicClusterCollection> bClusters1;
94			edm::Handle<reco::BasicClusterCollection> bClusters2;
95
96
97	yilmaz	1.1	typedef vector<EcalRecHit>::const_iterator EcalIterator;
98	yilmaz	1.2	typedef vector<HFRecHit>::const_iterator HFIterator;
99			typedef vector<HBHERecHit>::const_iterator HBHEIterator;
100	yilmaz	1.1
101			edm::Handle<reco::CaloJetCollection> signalJets;
102
103	yilmaz	1.2	edm::InputTag HcalRecHitHFSrc1_;
104			edm::InputTag HcalRecHitHFSrc2_;
105			edm::InputTag HcalRecHitHBHESrc1_;
106			edm::InputTag HcalRecHitHBHESrc2_;
107			edm::InputTag EBSrc1_;
108			edm::InputTag EBSrc2_;
109			edm::InputTag EESrc1_;
110			edm::InputTag EESrc2_;
111	yilmaz	1.7
112			edm::InputTag BCSrc1_;
113			edm::InputTag BCSrc2_;
114
115	yilmaz	1.3	edm::InputTag signalTag_;
116	yilmaz	1.2
117	yilmaz	1.7	TNtuple* ntBC;
118	yilmaz	1.2	TNtuple* ntEB;
119			TNtuple* ntEE;
120			TNtuple* ntHBHE;
121			TNtuple* ntHF;
122	yilmaz	1.1	TNtuple* ntjet;
123
124			double cone;
125	yilmaz	1.5	bool jetsOnly_;
126	yilmaz	1.7	bool doBasicClusters_;
127	yilmaz	1.9	bool doJetCone_;
128	yilmaz	1.1
129			edm::Service<TFileService> fs;
130	yilmaz	1.6	CentralityProvider* centrality_;
131	yilmaz	1.1	const CaloGeometry *geo;
132			};
133
134			//
135			// constants, enums and typedefs
136			//
137
138			//
139			// static data member definitions
140			//
141
142			//
143			// constructors and destructor
144			//
145			RecHitComparison::RecHitComparison(const edm::ParameterSet& iConfig) :
146	yilmaz	1.6	centrality_(0),
147	yilmaz	1.1	geo(0),
148			cone(0.5)
149			{
150			//now do what ever initialization is needed
151	yilmaz	1.5	jetsOnly_ = iConfig.getUntrackedParameter<bool>("jetsOnly",false);
152	yilmaz	1.7	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
153	yilmaz	1.9
154			doJetCone_ = iConfig.getUntrackedParameter<bool>("doJetCone",false);
155	yilmaz	1.3	signalTag_ = iConfig.getUntrackedParameter<edm::InputTag>("signalJets",edm::InputTag("iterativeCone5CaloJets","","SIGNAL"));
156
157	yilmaz	1.9	if(!doJetCone_) jetsOnly_ = 0;
158
159	yilmaz	1.2	HcalRecHitHFSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc1",edm::InputTag("hfreco"));
160			HcalRecHitHFSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc2",edm::InputTag("hfreco"));
161			HcalRecHitHBHESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc1",edm::InputTag("hbhereco"));
162			HcalRecHitHBHESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc2",edm::InputTag("hbhereco"));
163	yilmaz	1.10	EBSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECOBKG"));
164			EBSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","S"));
165	yilmaz	1.2	EESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEE","RECO"));
166			EESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEE","SIGNALONLY"));
167	yilmaz	1.7	BCSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
168			BCSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
169	yilmaz	1.9
170	yilmaz	1.1	}
171
172
173			RecHitComparison::~RecHitComparison()
174			{
175
176			// do anything here that needs to be done at desctruction time
177			// (e.g. close files, deallocate resources etc.)
178
179			}
180
181
182			//
183			// member functions
184			//
185
186			// ------------ method called to for each event ------------
187			void
188			RecHitComparison::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
189			{
190	yilmaz	1.6	if(!centrality_) centrality_ = new CentralityProvider(iSetup);
191	yilmaz	1.1	if(!geo){
192			edm::ESHandle<CaloGeometry> pGeo;
193			iSetup.get<CaloGeometryRecord>().get(pGeo);
194			geo = pGeo.product();
195			}
196
197			using namespace edm;
198	yilmaz	1.2	ev.getByLabel(EBSrc1_,ebHits1);
199	yilmaz	1.6	ev.getByLabel(EBSrc2_,ebHits2);
200	yilmaz	1.9
201			if(doJetCone_) ev.getByLabel(signalTag_,signalJets);
202	yilmaz	1.1
203	yilmaz	1.2	ev.getByLabel(HcalRecHitHFSrc1_,hfHits1);
204			ev.getByLabel(HcalRecHitHFSrc2_,hfHits2);
205			ev.getByLabel(HcalRecHitHBHESrc1_,hbheHits1);
206			ev.getByLabel(HcalRecHitHBHESrc2_,hbheHits2);
207			ev.getByLabel(EESrc1_,eeHits1);
208			ev.getByLabel(EESrc2_,eeHits2);
209
210	yilmaz	1.7	if(doBasicClusters_){
211			ev.getByLabel(BCSrc1_,bClusters1);
212			ev.getByLabel(BCSrc2_,bClusters2);
213			}
214
215	yilmaz	1.6	centrality_->newEvent(ev,iSetup);
216	yilmaz	1.10
217	yilmaz	1.8	double hf = centrality_->centralityValue();
218			int bin = centrality_->getBin();
219	yilmaz	1.1
220			vector<double> fFull;
221			vector<double> f05;
222			vector<double> f1;
223			vector<double> f15;
224			vector<double> f2;
225			vector<double> f25;
226			vector<double> f3;
227
228	yilmaz	1.10	int njets = 0;
229
230			if(doJetCone_) njets = signalJets->size();
231	yilmaz	1.1	fFull.reserve(njets);
232			f05.reserve(njets);
233			f1.reserve(njets);
234			f15.reserve(njets);
235			f2.reserve(njets);
236			f25.reserve(njets);
237			f3.reserve(njets);
238
239	yilmaz	1.10	if(doJetCone_){
240			for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
241			fFull.push_back(0);
242			f05.push_back(0);
243			f1.push_back(0);
244			f15.push_back(0);
245			f2.push_back(0);
246			f25.push_back(0);
247			f3.push_back(0);
248			}
249	yilmaz	1.1	}
250
251	yilmaz	1.2	for(unsigned int j1 = 0 ; j1 < ebHits1->size(); ++j1){
252	yilmaz	1.1
253	yilmaz	1.2	const EcalRecHit& jet1 = (*ebHits1)[j1];
254	yilmaz	1.1	double e1 = jet1.energy();
255
256			const GlobalPoint& pos1=geo->getPosition(jet1.id());
257			double eta1 = pos1.eta();
258			double phi1 = pos1.phi();
259			double et1 = e1*sin(pos1.theta());
260
261			double drjet = -1;
262			double jetpt = -1;
263			bool isjet = false;
264			int matchedJet = -1;
265
266	yilmaz	1.9	if(doJetCone_){
267			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
268			const reco::CaloJet & jet = (*signalJets)[j];
269			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
270			if(dr < cone){
271	yilmaz	1.1	jetpt = jet.pt();
272			drjet = dr;
273			isjet = true;
274			matchedJet = j;
275			fFull[j] += et1;
276	yilmaz	1.9
277	yilmaz	1.1	if(et1 > 0.5){
278	yilmaz	1.9	f05[j] += et1;
279	yilmaz	1.1	}
280			if(et1 > 1.){
281			f1[j] += et1;
282			}
283			if(et1 > 1.5){
284			f15[j] += et1;
285			}
286			if(et1 > 2){
287			f2[j] += et1;
288			}
289			if(et1 > 2.5){
290			f25[j] += et1;
291			}
292			if(et1 > 3.){
293			f3[j] += et1;
294			}
295	yilmaz	1.9	}
296			}
297			}
298	yilmaz	1.1
299			GlobalPoint pos2;
300			double e2 = -1;
301	yilmaz	1.2	EcalIterator hitit = ebHits2->find(jet1.id());
302			if(hitit != ebHits2->end()){
303	yilmaz	1.1	e2 = hitit->energy();
304			pos2=geo->getPosition(hitit->id());
305			}else{
306			e2 = 0;
307			pos2 = pos1;
308			}
309
310			double eta2 = pos2.eta();
311			double phi2 = pos2.eta();
312			double et2 = e2*sin(pos2.theta());
313	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntEB->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
314	yilmaz	1.2	}
315
316			for(unsigned int i = 0; i < eeHits1->size(); ++i){
317			const EcalRecHit & jet1= (*eeHits1)[i];
318			double e1 = jet1.energy();
319			const GlobalPoint& pos1=geo->getPosition(jet1.id());
320			double eta1 = pos1.eta();
321			double phi1 = pos1.phi();
322			double et1 = e1*sin(pos1.theta());
323			double drjet = -1;
324			double jetpt = -1;
325			bool isjet = false;
326			int matchedJet = -1;
327	yilmaz	1.9	if(doJetCone_){
328			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
329			const reco::CaloJet & jet = (*signalJets)[j];
330			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
331			if(dr < cone){
332			jetpt = jet.pt();
333			drjet = dr;
334			isjet = true;
335			matchedJet = j;
336			}
337	yilmaz	1.2	}
338			}
339	yilmaz	1.9
340	yilmaz	1.2	GlobalPoint pos2;
341			double e2 = -1;
342			EcalIterator hitit = eeHits2->find(jet1.id());
343			if(hitit != eeHits2->end()){
344			e2 = hitit->energy();
345			pos2=geo->getPosition(hitit->id());
346			}else{
347			e2 = 0;
348			pos2 = pos1;
349			}
350			double eta2 = pos2.eta();
351			double phi2 = pos2.eta();
352			double et2 = e2*sin(pos2.theta());
353	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntEE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
354	yilmaz	1.2	}
355
356			for(unsigned int i = 0; i < hbheHits1->size(); ++i){
357			const HBHERecHit & jet1= (*hbheHits1)[i];
358			double e1 = jet1.energy();
359			const GlobalPoint& pos1=geo->getPosition(jet1.id());
360			double eta1 = pos1.eta();
361			double phi1 = pos1.phi();
362			double et1 = e1*sin(pos1.theta());
363			double drjet = -1;
364			double jetpt = -1;
365			bool isjet = false;
366			int matchedJet = -1;
367	yilmaz	1.9	if(doJetCone_){
368			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
369			const reco::CaloJet & jet = (*signalJets)[j];
370			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
371			if(dr < cone){
372			jetpt = jet.pt();
373			drjet = dr;
374			isjet = true;
375			matchedJet = j;
376			}
377	yilmaz	1.2	}
378			}
379	yilmaz	1.9
380	yilmaz	1.2	GlobalPoint pos2;
381			double e2 = -1;
382			HBHEIterator hitit = hbheHits2->find(jet1.id());
383			if(hitit != hbheHits2->end()){
384			e2 = hitit->energy();
385			pos2=geo->getPosition(hitit->id());
386			}else{
387			e2 = 0;
388			pos2 = pos1;
389			}
390			double eta2 = pos2.eta();
391			double phi2 = pos2.eta();
392			double et2 = e2*sin(pos2.theta());
393	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntHBHE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
394	yilmaz	1.2	}
395
396			for(unsigned int i = 0; i < hfHits1->size(); ++i){
397			const HFRecHit & jet1= (*hfHits1)[i];
398			double e1 = jet1.energy();
399			const GlobalPoint& pos1=geo->getPosition(jet1.id());
400			double eta1 = pos1.eta();
401			double phi1 = pos1.phi();
402			double et1 = e1*sin(pos1.theta());
403			double drjet = -1;
404			double jetpt = -1;
405			bool isjet = false;
406			int matchedJet = -1;
407	yilmaz	1.9	if(doJetCone_){
408			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
409			const reco::CaloJet & jet = (*signalJets)[j];
410			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
411			if(dr < cone){
412			jetpt = jet.pt();
413			drjet = dr;
414			isjet = true;
415			matchedJet = j;
416			}
417	yilmaz	1.2	}
418			}
419			GlobalPoint pos2;
420			double e2 = -1;
421			HFIterator hitit = hfHits2->find(jet1.id());
422			if(hitit != hfHits2->end()){
423			e2 = hitit->energy();
424			pos2=geo->getPosition(hitit->id());
425			}else{
426			e2 = 0;
427			pos2 = pos1;
428			}
429			double eta2 = pos2.eta();
430			double phi2 = pos2.eta();
431			double et2 = e2*sin(pos2.theta());
432	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntHF->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
433	yilmaz	1.1	}
434
435	yilmaz	1.10	if(doJetCone_){
436			for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
437			const reco::CaloJet & jet = (*signalJets)[j1];
438			double em = (jet.emEnergyInEB() + jet.emEnergyInEE()) * sin(jet.theta());
439			double emf = jet.emEnergyFraction();
440			double pt = jet.pt();
441			double eta = jet.eta();
442			ntjet->Fill(bin,pt,eta,fFull[j1],f05[j1],f1[j1],f15[j1],f2[j1],f25[j1],f3[j1],em,emf);
443			}
444	yilmaz	1.1	}
445
446			}
447
448
449			// ------------ method called once each job just before starting event loop ------------
450			void
451			RecHitComparison::beginJob()
452			{
453	yilmaz	1.6	ntEB = fs->make<TNtuple>("ntEB","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
454			ntEE = fs->make<TNtuple>("ntEE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
455			ntHBHE = fs->make<TNtuple>("ntHBHE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
456			ntHF = fs->make<TNtuple>("ntHF","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
457	yilmaz	1.7
458			ntBC = fs->make<TNtuple>("ntBC","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
459
460	yilmaz	1.1	ntjet = fs->make<TNtuple>("ntjet","","bin:pt:eta:ethit:f05:f1:f15:f2:f25:f3:em:emf");
461
462			}
463
464			// ------------ method called once each job just after ending the event loop ------------
465			void
466			RecHitComparison::endJob() {
467			}
468
469			//define this as a plug-in
470			DEFINE_FWK_MODULE(RecHitComparison);