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.7 2010/11/01 21:48:31 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_;

   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);
   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"));
  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);
   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 = 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");

   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.8
Committed:	Tue Nov 2 13:19:06 2010 UTC (14 years, 6 months ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Changes since 1.7:	+3 -3 lines
Log Message:	for new centrality provider
#	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.8	// $Id: RecHitComparison.cc,v 1.7 2010/11/01 21:48:31 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.1
128			edm::Service<TFileService> fs;
129	yilmaz	1.6	CentralityProvider* centrality_;
130	yilmaz	1.1	const CaloGeometry *geo;
131			};
132
133			//
134			// constants, enums and typedefs
135			//
136
137			//
138			// static data member definitions
139			//
140
141			//
142			// constructors and destructor
143			//
144			RecHitComparison::RecHitComparison(const edm::ParameterSet& iConfig) :
145	yilmaz	1.6	centrality_(0),
146	yilmaz	1.1	geo(0),
147			cone(0.5)
148			{
149			//now do what ever initialization is needed
150	yilmaz	1.5	jetsOnly_ = iConfig.getUntrackedParameter<bool>("jetsOnly",false);
151	yilmaz	1.7	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
152	yilmaz	1.3	signalTag_ = iConfig.getUntrackedParameter<edm::InputTag>("signalJets",edm::InputTag("iterativeCone5CaloJets","","SIGNAL"));
153
154	yilmaz	1.2	HcalRecHitHFSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc1",edm::InputTag("hfreco"));
155			HcalRecHitHFSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc2",edm::InputTag("hfreco"));
156			HcalRecHitHBHESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc1",edm::InputTag("hbhereco"));
157			HcalRecHitHBHESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc2",edm::InputTag("hbhereco"));
158			EBSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
159			EBSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
160			EESrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEE","RECO"));
161			EESrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEE","SIGNALONLY"));
162	yilmaz	1.7	BCSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
163			BCSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
164	yilmaz	1.1	}
165
166
167			RecHitComparison::~RecHitComparison()
168			{
169
170			// do anything here that needs to be done at desctruction time
171			// (e.g. close files, deallocate resources etc.)
172
173			}
174
175
176			//
177			// member functions
178			//
179
180			// ------------ method called to for each event ------------
181			void
182			RecHitComparison::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
183			{
184	yilmaz	1.6	if(!centrality_) centrality_ = new CentralityProvider(iSetup);
185	yilmaz	1.1	if(!geo){
186			edm::ESHandle<CaloGeometry> pGeo;
187			iSetup.get<CaloGeometryRecord>().get(pGeo);
188			geo = pGeo.product();
189			}
190
191			using namespace edm;
192	yilmaz	1.2	ev.getByLabel(EBSrc1_,ebHits1);
193	yilmaz	1.6	ev.getByLabel(EBSrc2_,ebHits2);
194	yilmaz	1.3	ev.getByLabel(signalTag_,signalJets);
195	yilmaz	1.1
196	yilmaz	1.2	ev.getByLabel(HcalRecHitHFSrc1_,hfHits1);
197			ev.getByLabel(HcalRecHitHFSrc2_,hfHits2);
198			ev.getByLabel(HcalRecHitHBHESrc1_,hbheHits1);
199			ev.getByLabel(HcalRecHitHBHESrc2_,hbheHits2);
200
201			ev.getByLabel(EESrc1_,eeHits1);
202			ev.getByLabel(EESrc2_,eeHits2);
203
204	yilmaz	1.7	if(doBasicClusters_){
205			ev.getByLabel(BCSrc1_,bClusters1);
206			ev.getByLabel(BCSrc2_,bClusters2);
207			}
208
209	yilmaz	1.6	centrality_->newEvent(ev,iSetup);
210	yilmaz	1.8	double hf = centrality_->centralityValue();
211			int bin = centrality_->getBin();
212	yilmaz	1.1
213			vector<double> fFull;
214			vector<double> f05;
215			vector<double> f1;
216			vector<double> f15;
217			vector<double> f2;
218			vector<double> f25;
219			vector<double> f3;
220
221			int njets = signalJets->size();
222			fFull.reserve(njets);
223			f05.reserve(njets);
224			f1.reserve(njets);
225			f15.reserve(njets);
226			f2.reserve(njets);
227			f25.reserve(njets);
228			f3.reserve(njets);
229
230			for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
231			fFull.push_back(0);
232			f05.push_back(0);
233			f1.push_back(0);
234			f15.push_back(0);
235			f2.push_back(0);
236			f25.push_back(0);
237			f3.push_back(0);
238			}
239
240	yilmaz	1.2	for(unsigned int j1 = 0 ; j1 < ebHits1->size(); ++j1){
241	yilmaz	1.1
242	yilmaz	1.2	const EcalRecHit& jet1 = (*ebHits1)[j1];
243	yilmaz	1.1	double e1 = jet1.energy();
244
245			const GlobalPoint& pos1=geo->getPosition(jet1.id());
246			double eta1 = pos1.eta();
247			double phi1 = pos1.phi();
248			double et1 = e1*sin(pos1.theta());
249
250			double drjet = -1;
251			double jetpt = -1;
252			bool isjet = false;
253			int matchedJet = -1;
254
255			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
256			const reco::CaloJet & jet = (*signalJets)[j];
257			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
258			if(dr < cone){
259			jetpt = jet.pt();
260			drjet = dr;
261			isjet = true;
262			matchedJet = j;
263			fFull[j] += et1;
264
265			if(et1 > 0.5){
266			f05[j] += et1;
267			}
268			if(et1 > 1.){
269			f1[j] += et1;
270			}
271			if(et1 > 1.5){
272			f15[j] += et1;
273			}
274			if(et1 > 2){
275			f2[j] += et1;
276			}
277			if(et1 > 2.5){
278			f25[j] += et1;
279			}
280			if(et1 > 3.){
281			f3[j] += et1;
282			}
283			}
284			}
285
286			GlobalPoint pos2;
287			double e2 = -1;
288	yilmaz	1.2	EcalIterator hitit = ebHits2->find(jet1.id());
289			if(hitit != ebHits2->end()){
290	yilmaz	1.1	e2 = hitit->energy();
291			pos2=geo->getPosition(hitit->id());
292			}else{
293			e2 = 0;
294			pos2 = pos1;
295			}
296
297			double eta2 = pos2.eta();
298			double phi2 = pos2.eta();
299			double et2 = e2*sin(pos2.theta());
300	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntEB->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
301	yilmaz	1.2	}
302
303			for(unsigned int i = 0; i < eeHits1->size(); ++i){
304			const EcalRecHit & jet1= (*eeHits1)[i];
305			double e1 = jet1.energy();
306			const GlobalPoint& pos1=geo->getPosition(jet1.id());
307			double eta1 = pos1.eta();
308			double phi1 = pos1.phi();
309			double et1 = e1*sin(pos1.theta());
310			double drjet = -1;
311			double jetpt = -1;
312			bool isjet = false;
313			int matchedJet = -1;
314			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
315			const reco::CaloJet & jet = (*signalJets)[j];
316			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
317			if(dr < cone){
318			jetpt = jet.pt();
319			drjet = dr;
320			isjet = true;
321			matchedJet = j;
322			}
323			}
324			GlobalPoint pos2;
325			double e2 = -1;
326			EcalIterator hitit = eeHits2->find(jet1.id());
327			if(hitit != eeHits2->end()){
328			e2 = hitit->energy();
329			pos2=geo->getPosition(hitit->id());
330			}else{
331			e2 = 0;
332			pos2 = pos1;
333			}
334			double eta2 = pos2.eta();
335			double phi2 = pos2.eta();
336			double et2 = e2*sin(pos2.theta());
337	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntEE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
338	yilmaz	1.2	}
339
340			for(unsigned int i = 0; i < hbheHits1->size(); ++i){
341			const HBHERecHit & jet1= (*hbheHits1)[i];
342			double e1 = jet1.energy();
343			const GlobalPoint& pos1=geo->getPosition(jet1.id());
344			double eta1 = pos1.eta();
345			double phi1 = pos1.phi();
346			double et1 = e1*sin(pos1.theta());
347			double drjet = -1;
348			double jetpt = -1;
349			bool isjet = false;
350			int matchedJet = -1;
351			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
352			const reco::CaloJet & jet = (*signalJets)[j];
353			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
354			if(dr < cone){
355			jetpt = jet.pt();
356			drjet = dr;
357			isjet = true;
358			matchedJet = j;
359			}
360			}
361			GlobalPoint pos2;
362			double e2 = -1;
363			HBHEIterator hitit = hbheHits2->find(jet1.id());
364			if(hitit != hbheHits2->end()){
365			e2 = hitit->energy();
366			pos2=geo->getPosition(hitit->id());
367			}else{
368			e2 = 0;
369			pos2 = pos1;
370			}
371			double eta2 = pos2.eta();
372			double phi2 = pos2.eta();
373			double et2 = e2*sin(pos2.theta());
374	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntHBHE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
375	yilmaz	1.2	}
376
377			for(unsigned int i = 0; i < hfHits1->size(); ++i){
378			const HFRecHit & jet1= (*hfHits1)[i];
379			double e1 = jet1.energy();
380			const GlobalPoint& pos1=geo->getPosition(jet1.id());
381			double eta1 = pos1.eta();
382			double phi1 = pos1.phi();
383			double et1 = e1*sin(pos1.theta());
384			double drjet = -1;
385			double jetpt = -1;
386			bool isjet = false;
387			int matchedJet = -1;
388			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
389			const reco::CaloJet & jet = (*signalJets)[j];
390			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
391			if(dr < cone){
392			jetpt = jet.pt();
393			drjet = dr;
394			isjet = true;
395			matchedJet = j;
396			}
397			}
398			GlobalPoint pos2;
399			double e2 = -1;
400			HFIterator hitit = hfHits2->find(jet1.id());
401			if(hitit != hfHits2->end()){
402			e2 = hitit->energy();
403			pos2=geo->getPosition(hitit->id());
404			}else{
405			e2 = 0;
406			pos2 = pos1;
407			}
408			double eta2 = pos2.eta();
409			double phi2 = pos2.eta();
410			double et2 = e2*sin(pos2.theta());
411	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntHF->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
412	yilmaz	1.1	}
413
414			for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
415			const reco::CaloJet & jet = (*signalJets)[j1];
416			double em = (jet.emEnergyInEB() + jet.emEnergyInEE()) * sin(jet.theta());
417			double emf = jet.emEnergyFraction();
418			double pt = jet.pt();
419			double eta = jet.eta();
420			ntjet->Fill(bin,pt,eta,fFull[j1],f05[j1],f1[j1],f15[j1],f2[j1],f25[j1],f3[j1],em,emf);
421			}
422
423			}
424
425
426			// ------------ method called once each job just before starting event loop ------------
427			void
428			RecHitComparison::beginJob()
429			{
430	yilmaz	1.6	ntEB = fs->make<TNtuple>("ntEB","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
431			ntEE = fs->make<TNtuple>("ntEE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
432			ntHBHE = fs->make<TNtuple>("ntHBHE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
433			ntHF = fs->make<TNtuple>("ntHF","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
434	yilmaz	1.7
435			ntBC = fs->make<TNtuple>("ntBC","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
436
437	yilmaz	1.1	ntjet = fs->make<TNtuple>("ntjet","","bin:pt:eta:ethit:f05:f1:f15:f2:f25:f3:em:emf");
438
439			}
440
441			// ------------ method called once each job just after ending the event loop ------------
442			void
443			RecHitComparison::endJob() {
444			}
445
446			//define this as a plug-in
447			DEFINE_FWK_MODULE(RecHitComparison);