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.4 2010/10/20 15:51:18 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/Centrality.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<reco::Centrality> cent;
   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_;
   edm::InputTag centTag_;

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

   double cone;
   bool jetsOnly_;

   edm::Service<TFileService> fs;
   const CentralityBins * cbins_;
   const CaloGeometry *geo;
};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
RecHitComparison::RecHitComparison(const edm::ParameterSet& iConfig) :
   cbins_(0),
   geo(0),
   cone(0.5)
{
   //now do what ever initialization is needed
   centTag_ =  iConfig.getUntrackedParameter<edm::InputTag>("centrality",edm::InputTag("hiCentrality","","RECO"));

   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(!cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
   if(!geo){
      edm::ESHandle<CaloGeometry> pGeo;
      iSetup.get<CaloGeometryRecord>().get(pGeo);
      geo = pGeo.product();
   }

   using namespace edm;

   ev.getByLabel(centTag_,cent);
   ev.getByLabel(EBSrc1_,ebHits1);
   ev.getByLabel(EBSrc1_,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);

   double hf = cent->EtHFhitSum();
   double sumet = cent->EtMidRapiditySum();
   int run = ev.id().run();
   run = 1;
   int bin = cbins_->getBin(hf);
   int margin = 0;

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