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.8 2010/11/02 13:19:06 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","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);

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

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

   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.9
Committed:	Wed Mar 30 12:13:22 2011 UTC (14 years, 1 month ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Changes since 1.8:	+51 -34 lines
Log Message:	BF
#	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.9	// $Id: RecHitComparison.cc,v 1.8 2010/11/02 13:19:06 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			EBSrc1_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
164			EBSrc2_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc2",edm::InputTag("ecalRecHit","EcalRecHitsEB","SIGNALONLY"));
165			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
208			ev.getByLabel(EESrc1_,eeHits1);
209			ev.getByLabel(EESrc2_,eeHits2);
210
211	yilmaz	1.7	if(doBasicClusters_){
212			ev.getByLabel(BCSrc1_,bClusters1);
213			ev.getByLabel(BCSrc2_,bClusters2);
214			}
215
216	yilmaz	1.6	centrality_->newEvent(ev,iSetup);
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			int njets = signalJets->size();
229			fFull.reserve(njets);
230			f05.reserve(njets);
231			f1.reserve(njets);
232			f15.reserve(njets);
233			f2.reserve(njets);
234			f25.reserve(njets);
235			f3.reserve(njets);
236
237			for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
238			fFull.push_back(0);
239			f05.push_back(0);
240			f1.push_back(0);
241			f15.push_back(0);
242			f2.push_back(0);
243			f25.push_back(0);
244			f3.push_back(0);
245			}
246
247	yilmaz	1.2	for(unsigned int j1 = 0 ; j1 < ebHits1->size(); ++j1){
248	yilmaz	1.1
249	yilmaz	1.2	const EcalRecHit& jet1 = (*ebHits1)[j1];
250	yilmaz	1.1	double e1 = jet1.energy();
251
252			const GlobalPoint& pos1=geo->getPosition(jet1.id());
253			double eta1 = pos1.eta();
254			double phi1 = pos1.phi();
255			double et1 = e1*sin(pos1.theta());
256
257			double drjet = -1;
258			double jetpt = -1;
259			bool isjet = false;
260			int matchedJet = -1;
261
262	yilmaz	1.9	if(doJetCone_){
263			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
264			const reco::CaloJet & jet = (*signalJets)[j];
265			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
266			if(dr < cone){
267	yilmaz	1.1	jetpt = jet.pt();
268			drjet = dr;
269			isjet = true;
270			matchedJet = j;
271			fFull[j] += et1;
272	yilmaz	1.9
273	yilmaz	1.1	if(et1 > 0.5){
274	yilmaz	1.9	f05[j] += et1;
275	yilmaz	1.1	}
276			if(et1 > 1.){
277			f1[j] += et1;
278			}
279			if(et1 > 1.5){
280			f15[j] += et1;
281			}
282			if(et1 > 2){
283			f2[j] += et1;
284			}
285			if(et1 > 2.5){
286			f25[j] += et1;
287			}
288			if(et1 > 3.){
289			f3[j] += et1;
290			}
291	yilmaz	1.9	}
292			}
293			}
294	yilmaz	1.1
295			GlobalPoint pos2;
296			double e2 = -1;
297	yilmaz	1.2	EcalIterator hitit = ebHits2->find(jet1.id());
298			if(hitit != ebHits2->end()){
299	yilmaz	1.1	e2 = hitit->energy();
300			pos2=geo->getPosition(hitit->id());
301			}else{
302			e2 = 0;
303			pos2 = pos1;
304			}
305
306			double eta2 = pos2.eta();
307			double phi2 = pos2.eta();
308			double et2 = e2*sin(pos2.theta());
309	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntEB->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
310	yilmaz	1.2	}
311
312			for(unsigned int i = 0; i < eeHits1->size(); ++i){
313			const EcalRecHit & jet1= (*eeHits1)[i];
314			double e1 = jet1.energy();
315			const GlobalPoint& pos1=geo->getPosition(jet1.id());
316			double eta1 = pos1.eta();
317			double phi1 = pos1.phi();
318			double et1 = e1*sin(pos1.theta());
319			double drjet = -1;
320			double jetpt = -1;
321			bool isjet = false;
322			int matchedJet = -1;
323	yilmaz	1.9	if(doJetCone_){
324			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
325			const reco::CaloJet & jet = (*signalJets)[j];
326			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
327			if(dr < cone){
328			jetpt = jet.pt();
329			drjet = dr;
330			isjet = true;
331			matchedJet = j;
332			}
333	yilmaz	1.2	}
334			}
335	yilmaz	1.9
336	yilmaz	1.2	GlobalPoint pos2;
337			double e2 = -1;
338			EcalIterator hitit = eeHits2->find(jet1.id());
339			if(hitit != eeHits2->end()){
340			e2 = hitit->energy();
341			pos2=geo->getPosition(hitit->id());
342			}else{
343			e2 = 0;
344			pos2 = pos1;
345			}
346			double eta2 = pos2.eta();
347			double phi2 = pos2.eta();
348			double et2 = e2*sin(pos2.theta());
349	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntEE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
350	yilmaz	1.2	}
351
352			for(unsigned int i = 0; i < hbheHits1->size(); ++i){
353			const HBHERecHit & jet1= (*hbheHits1)[i];
354			double e1 = jet1.energy();
355			const GlobalPoint& pos1=geo->getPosition(jet1.id());
356			double eta1 = pos1.eta();
357			double phi1 = pos1.phi();
358			double et1 = e1*sin(pos1.theta());
359			double drjet = -1;
360			double jetpt = -1;
361			bool isjet = false;
362			int matchedJet = -1;
363	yilmaz	1.9	if(doJetCone_){
364			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
365			const reco::CaloJet & jet = (*signalJets)[j];
366			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
367			if(dr < cone){
368			jetpt = jet.pt();
369			drjet = dr;
370			isjet = true;
371			matchedJet = j;
372			}
373	yilmaz	1.2	}
374			}
375	yilmaz	1.9
376	yilmaz	1.2	GlobalPoint pos2;
377			double e2 = -1;
378			HBHEIterator hitit = hbheHits2->find(jet1.id());
379			if(hitit != hbheHits2->end()){
380			e2 = hitit->energy();
381			pos2=geo->getPosition(hitit->id());
382			}else{
383			e2 = 0;
384			pos2 = pos1;
385			}
386			double eta2 = pos2.eta();
387			double phi2 = pos2.eta();
388			double et2 = e2*sin(pos2.theta());
389	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntHBHE->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
390	yilmaz	1.2	}
391
392			for(unsigned int i = 0; i < hfHits1->size(); ++i){
393			const HFRecHit & jet1= (*hfHits1)[i];
394			double e1 = jet1.energy();
395			const GlobalPoint& pos1=geo->getPosition(jet1.id());
396			double eta1 = pos1.eta();
397			double phi1 = pos1.phi();
398			double et1 = e1*sin(pos1.theta());
399			double drjet = -1;
400			double jetpt = -1;
401			bool isjet = false;
402			int matchedJet = -1;
403	yilmaz	1.9	if(doJetCone_){
404			for(unsigned int j = 0 ; j < signalJets->size(); ++j){
405			const reco::CaloJet & jet = (*signalJets)[j];
406			double dr = reco::deltaR(eta1,phi1,jet.eta(),jet.phi());
407			if(dr < cone){
408			jetpt = jet.pt();
409			drjet = dr;
410			isjet = true;
411			matchedJet = j;
412			}
413	yilmaz	1.2	}
414			}
415			GlobalPoint pos2;
416			double e2 = -1;
417			HFIterator hitit = hfHits2->find(jet1.id());
418			if(hitit != hfHits2->end()){
419			e2 = hitit->energy();
420			pos2=geo->getPosition(hitit->id());
421			}else{
422			e2 = 0;
423			pos2 = pos1;
424			}
425			double eta2 = pos2.eta();
426			double phi2 = pos2.eta();
427			double et2 = e2*sin(pos2.theta());
428	yilmaz	1.6	if(!jetsOnly_ \|\| isjet) ntHF->Fill(e1,et1,e2,et2,eta2,phi2,hf,bin,jetpt,drjet);
429	yilmaz	1.1	}
430
431			for(unsigned int j1 = 0 ; j1 < signalJets->size(); ++j1){
432			const reco::CaloJet & jet = (*signalJets)[j1];
433			double em = (jet.emEnergyInEB() + jet.emEnergyInEE()) * sin(jet.theta());
434			double emf = jet.emEnergyFraction();
435			double pt = jet.pt();
436			double eta = jet.eta();
437			ntjet->Fill(bin,pt,eta,fFull[j1],f05[j1],f1[j1],f15[j1],f2[j1],f25[j1],f3[j1],em,emf);
438			}
439
440			}
441
442
443			// ------------ method called once each job just before starting event loop ------------
444			void
445			RecHitComparison::beginJob()
446			{
447	yilmaz	1.6	ntEB = fs->make<TNtuple>("ntEB","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
448			ntEE = fs->make<TNtuple>("ntEE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
449			ntHBHE = fs->make<TNtuple>("ntHBHE","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
450			ntHF = fs->make<TNtuple>("ntHF","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
451	yilmaz	1.7
452			ntBC = fs->make<TNtuple>("ntBC","","e1:et1:e2:et2:eta:phi:hf:bin:ptjet:drjet");
453
454	yilmaz	1.1	ntjet = fs->make<TNtuple>("ntjet","","bin:pt:eta:ethit:f05:f1:f15:f2:f25:f3:em:emf");
455
456			}
457
458			// ------------ method called once each job just after ending the event loop ------------
459			void
460			RecHitComparison::endJob() {
461			}
462
463			//define this as a plug-in
464			DEFINE_FWK_MODULE(RecHitComparison);