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.5 2010/10/22 13:34:09 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 "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;

   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 signalTag_;

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

   double cone;
   bool jetsOnly_;

   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);
   signalTag_ = iConfig.getUntrackedParameter<edm::InputTag>("signalJets",edm::InputTag("iterativeCone5CaloJets","","SIGNAL"));

  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","RECO"));
  EBSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
  EESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEE","RECO"));
  EESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEE","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);
   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);

   centrality_->newEvent(ev,iSetup);
   double hf = centrality_->centralityValue(ev);
   int bin = centrality_->getBin(ev);

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

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

   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;

      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;
     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;
     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;
     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);
   }

   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");
   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.6
Committed:	Sun Oct 24 14:39:32 2010 UTC (14 years, 6 months ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Changes since 1.5:	+17 -29 lines
Log Message:	Some fixes
#	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.6	// $Id: RecHitComparison.cc,v 1.5 2010/10/22 13:34:09 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			#include "TNtuple.h"
59
60			using namespace std;
61
62			//
63			// class declaration
64			//
65
66			class RecHitComparison : public edm::EDAnalyzer {
67			public:
68			explicit RecHitComparison(const edm::ParameterSet&);
69			~RecHitComparison();
70
71
72			private:
73			virtual void beginJob() ;
74			virtual void analyze(const edm::Event&, const edm::EventSetup&);
75			virtual void endJob() ;
76
77			// ----------member data ---------------------------
78			edm::Handle<vector<double> > ktRhos;
79			edm::Handle<vector<double> > akRhos;
80
81	yilmaz	1.2	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits1;
82			edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits2;
83			edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits1;
84			edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits2;
85
86			edm::Handle<HFRecHitCollection> hfHits1;
87			edm::Handle<HFRecHitCollection> hfHits2;
88			edm::Handle<HBHERecHitCollection> hbheHits1;
89			edm::Handle<HBHERecHitCollection> hbheHits2;
90	yilmaz	1.1
91			typedef vector<EcalRecHit>::const_iterator EcalIterator;
92	yilmaz	1.2	typedef vector<HFRecHit>::const_iterator HFIterator;
93			typedef vector<HBHERecHit>::const_iterator HBHEIterator;
94	yilmaz	1.1
95			edm::Handle<reco::CaloJetCollection> signalJets;
96
97	yilmaz	1.2	edm::InputTag HcalRecHitHFSrc1_;
98			edm::InputTag HcalRecHitHFSrc2_;
99			edm::InputTag HcalRecHitHBHESrc1_;
100			edm::InputTag HcalRecHitHBHESrc2_;
101			edm::InputTag EBSrc1_;
102			edm::InputTag EBSrc2_;
103			edm::InputTag EESrc1_;
104			edm::InputTag EESrc2_;
105	yilmaz	1.3	edm::InputTag signalTag_;
106	yilmaz	1.2
107			TNtuple* ntEB;
108			TNtuple* ntEE;
109			TNtuple* ntHBHE;
110			TNtuple* ntHF;
111	yilmaz	1.1	TNtuple* ntjet;
112
113			double cone;
114	yilmaz	1.5	bool jetsOnly_;
115	yilmaz	1.1
116			edm::Service<TFileService> fs;
117	yilmaz	1.6	CentralityProvider* centrality_;
118	yilmaz	1.1	const CaloGeometry *geo;
119			};
120
121			//
122			// constants, enums and typedefs
123			//
124
125			//
126			// static data member definitions
127			//
128
129			//
130			// constructors and destructor
131			//
132			RecHitComparison::RecHitComparison(const edm::ParameterSet& iConfig) :
133	yilmaz	1.6	centrality_(0),
134	yilmaz	1.1	geo(0),
135			cone(0.5)
136			{
137			//now do what ever initialization is needed
138	yilmaz	1.5	jetsOnly_ = iConfig.getUntrackedParameter<bool>("jetsOnly",false);
139	yilmaz	1.3	signalTag_ = iConfig.getUntrackedParameter<edm::InputTag>("signalJets",edm::InputTag("iterativeCone5CaloJets","","SIGNAL"));
140
141	yilmaz	1.2	HcalRecHitHFSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc1",edm::InputTag("hfreco"));
142			HcalRecHitHFSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc2",edm::InputTag("hfreco"));
143			HcalRecHitHBHESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc1",edm::InputTag("hbhereco"));
144			HcalRecHitHBHESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc2",edm::InputTag("hbhereco"));
145			EBSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
146			EBSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
147			EESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEE","RECO"));
148			EESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEE","SIGNALONLY"));
149	yilmaz	1.1
150			}
151
152
153			RecHitComparison::~RecHitComparison()
154			{
155
156			// do anything here that needs to be done at desctruction time
157			// (e.g. close files, deallocate resources etc.)
158
159			}
160
161
162			//
163			// member functions
164			//
165
166			// ------------ method called to for each event ------------
167			void
168			RecHitComparison::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
169			{
170	yilmaz	1.6	if(!centrality_) centrality_ = new CentralityProvider(iSetup);
171	yilmaz	1.1	if(!geo){
172			edm::ESHandle<CaloGeometry> pGeo;
173			iSetup.get<CaloGeometryRecord>().get(pGeo);
174			geo = pGeo.product();
175			}
176
177			using namespace edm;
178	yilmaz	1.2	ev.getByLabel(EBSrc1_,ebHits1);
179	yilmaz	1.6	ev.getByLabel(EBSrc2_,ebHits2);
180	yilmaz	1.3	ev.getByLabel(signalTag_,signalJets);
181	yilmaz	1.1
182	yilmaz	1.2	ev.getByLabel(HcalRecHitHFSrc1_,hfHits1);
183			ev.getByLabel(HcalRecHitHFSrc2_,hfHits2);
184			ev.getByLabel(HcalRecHitHBHESrc1_,hbheHits1);
185			ev.getByLabel(HcalRecHitHBHESrc2_,hbheHits2);
186
187			ev.getByLabel(EESrc1_,eeHits1);
188			ev.getByLabel(EESrc2_,eeHits2);
189
190	yilmaz	1.6	centrality_->newEvent(ev,iSetup);
191			double hf = centrality_->centralityValue(ev);
192			int bin = centrality_->getBin(ev);
193	yilmaz	1.1
194			vector<double> fFull;
195			vector<double> f05;
196			vector<double> f1;
197			vector<double> f15;
198			vector<double> f2;
199			vector<double> f25;
200			vector<double> f3;
201
202			int njets = signalJets->size();
203			fFull.reserve(njets);
204			f05.reserve(njets);
205			f1.reserve(njets);
206			f15.reserve(njets);
207			f2.reserve(njets);
208			f25.reserve(njets);
209			f3.reserve(njets);
210
211			for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
212			fFull.push_back(0);
213			f05.push_back(0);
214			f1.push_back(0);
215			f15.push_back(0);
216			f2.push_back(0);
217			f25.push_back(0);
218			f3.push_back(0);
219			}
220
221	yilmaz	1.2	for(unsigned int j1 = 0 ; j1 < ebHits1->size(); ++j1){
222	yilmaz	1.1
223	yilmaz	1.2	const EcalRecHit& jet1 = (*ebHits1)[j1];
224	yilmaz	1.1	double e1 = jet1.energy();
225
226			const GlobalPoint& pos1=geo->getPosition(jet1.id());
227			double eta1 = pos1.eta();
228			double phi1 = pos1.phi();
229			double et1 = e1*sin(pos1.theta());
230
231			double drjet = -1;
232			double jetpt = -1;
233			bool isjet = false;
234			int matchedJet = -1;
235
236			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
237			const reco::CaloJet & jet = (*signalJets)[j];
238			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
239			if(dr < cone){
240			jetpt = jet.pt();
241			drjet = dr;
242			isjet = true;
243			matchedJet = j;
244			fFull[j] += et1;
245
246			if(et1 > 0.5){
247			f05[j] += et1;
248			}
249			if(et1 > 1.){
250			f1[j] += et1;
251			}
252			if(et1 > 1.5){
253			f15[j] += et1;
254			}
255			if(et1 > 2){
256			f2[j] += et1;
257			}
258			if(et1 > 2.5){
259			f25[j] += et1;
260			}
261			if(et1 > 3.){
262			f3[j] += et1;
263			}
264			}
265			}
266
267			GlobalPoint pos2;
268			double e2 = -1;
269	yilmaz	1.2	EcalIterator hitit = ebHits2->find(jet1.id());
270			if(hitit != ebHits2->end()){
271	yilmaz	1.1	e2 = hitit->energy();
272			pos2=geo->getPosition(hitit->id());
273			}else{
274			e2 = 0;
275			pos2 = pos1;
276			}
277
278			double eta2 = pos2.eta();
279			double phi2 = pos2.eta();
280			double et2 = e2*sin(pos2.theta());
281	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntEB->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
282	yilmaz	1.2	}
283
284			for(unsigned int i = 0; i < eeHits1->size(); ++i){
285			const EcalRecHit & jet1= (*eeHits1)[i];
286			double e1 = jet1.energy();
287			const GlobalPoint& pos1=geo->getPosition(jet1.id());
288			double eta1 = pos1.eta();
289			double phi1 = pos1.phi();
290			double et1 = e1*sin(pos1.theta());
291			double drjet = -1;
292			double jetpt = -1;
293			bool isjet = false;
294			int matchedJet = -1;
295			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
296			const reco::CaloJet & jet = (*signalJets)[j];
297			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
298			if(dr < cone){
299			jetpt = jet.pt();
300			drjet = dr;
301			isjet = true;
302			matchedJet = j;
303			}
304			}
305			GlobalPoint pos2;
306			double e2 = -1;
307			EcalIterator hitit = eeHits2->find(jet1.id());
308			if(hitit != eeHits2->end()){
309			e2 = hitit->energy();
310			pos2=geo->getPosition(hitit->id());
311			}else{
312			e2 = 0;
313			pos2 = pos1;
314			}
315			double eta2 = pos2.eta();
316			double phi2 = pos2.eta();
317			double et2 = e2*sin(pos2.theta());
318	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntEE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
319	yilmaz	1.2	}
320
321			for(unsigned int i = 0; i < hbheHits1->size(); ++i){
322			const HBHERecHit & jet1= (*hbheHits1)[i];
323			double e1 = jet1.energy();
324			const GlobalPoint& pos1=geo->getPosition(jet1.id());
325			double eta1 = pos1.eta();
326			double phi1 = pos1.phi();
327			double et1 = e1*sin(pos1.theta());
328			double drjet = -1;
329			double jetpt = -1;
330			bool isjet = false;
331			int matchedJet = -1;
332			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
333			const reco::CaloJet & jet = (*signalJets)[j];
334			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
335			if(dr < cone){
336			jetpt = jet.pt();
337			drjet = dr;
338			isjet = true;
339			matchedJet = j;
340			}
341			}
342			GlobalPoint pos2;
343			double e2 = -1;
344			HBHEIterator hitit = hbheHits2->find(jet1.id());
345			if(hitit != hbheHits2->end()){
346			e2 = hitit->energy();
347			pos2=geo->getPosition(hitit->id());
348			}else{
349			e2 = 0;
350			pos2 = pos1;
351			}
352			double eta2 = pos2.eta();
353			double phi2 = pos2.eta();
354			double et2 = e2*sin(pos2.theta());
355	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntHBHE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
356	yilmaz	1.2	}
357
358			for(unsigned int i = 0; i < hfHits1->size(); ++i){
359			const HFRecHit & jet1= (*hfHits1)[i];
360			double e1 = jet1.energy();
361			const GlobalPoint& pos1=geo->getPosition(jet1.id());
362			double eta1 = pos1.eta();
363			double phi1 = pos1.phi();
364			double et1 = e1*sin(pos1.theta());
365			double drjet = -1;
366			double jetpt = -1;
367			bool isjet = false;
368			int matchedJet = -1;
369			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
370			const reco::CaloJet & jet = (*signalJets)[j];
371			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
372			if(dr < cone){
373			jetpt = jet.pt();
374			drjet = dr;
375			isjet = true;
376			matchedJet = j;
377			}
378			}
379			GlobalPoint pos2;
380			double e2 = -1;
381			HFIterator hitit = hfHits2->find(jet1.id());
382			if(hitit != hfHits2->end()){
383			e2 = hitit->energy();
384			pos2=geo->getPosition(hitit->id());
385			}else{
386			e2 = 0;
387			pos2 = pos1;
388			}
389			double eta2 = pos2.eta();
390			double phi2 = pos2.eta();
391			double et2 = e2*sin(pos2.theta());
392	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntHF->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
393	yilmaz	1.1	}
394
395			for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
396			const reco::CaloJet & jet = (*signalJets)[j1];
397			double em = (jet.emEnergyInEB() + jet.emEnergyInEE()) * sin(jet.theta());
398			double emf = jet.emEnergyFraction();
399			double pt = jet.pt();
400			double eta = jet.eta();
401			ntjet->Fill(bin,pt,eta,fFull[j1],f05[j1],f1[j1],f15[j1],f2[j1],f25[j1],f3[j1],em,emf);
402			}
403
404			}
405
406
407			// ------------ method called once each job just before starting event loop ------------
408			void
409			RecHitComparison::beginJob()
410			{
411	yilmaz	1.6	ntEB = fs->make<TNtuple>("ntEB","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
412			ntEE = fs->make<TNtuple>("ntEE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
413			ntHBHE = fs->make<TNtuple>("ntHBHE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
414			ntHF = fs->make<TNtuple>("ntHF","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
415	yilmaz	1.1	ntjet = fs->make<TNtuple>("ntjet","","bin:pt:eta:ethit:f05:f1:f15:f2:f25:f3:em:emf");
416
417			}
418
419			// ------------ method called once each job just after ending the event loop ------------
420			void
421			RecHitComparison::endJob() {
422			}
423
424			//define this as a plug-in
425			DEFINE_FWK_MODULE(RecHitComparison);