JetAnalysis/src/RecHitTreeProducer.cc

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

 Description: [one line class summary]

 Implementation:
     [Notes on implementation]
*/
//
// Original Author:  Yetkin Yilmaz
//         Created:  Tue Sep  7 11:38:19 EDT 2010
// $Id: RecHitTreeProducer.cc,v 1.3 2010/10/20 15:01:11 nart 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/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
#include "DataFormats/DetId/interface/DetId.h"


#include "TNtuple.h"

using namespace std;

#define MAXHITS 1000000

struct MyRecHit{

  int n;

  float e[MAXHITS];
  float et[MAXHITS];
  float eta[MAXHITS];
  float phi[MAXHITS];
  bool isjet[MAXHITS];

};


//
// class declaration
//

class RecHitTreeProducer : public edm::EDAnalyzer {
   public:
      explicit RecHitTreeProducer(const edm::ParameterSet&);
      ~RecHitTreeProducer();
  double       getEt(const DetId &id, double energy);
  double       getEta(const DetId &id);
  double       getPhi(const DetId &id);


   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> > > ebHits;
  edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits;

  edm::Handle<HFRecHitCollection> hfHits;
  edm::Handle<HBHERecHitCollection> hbheHits;

  typedef vector<EcalRecHit>::const_iterator EcalIterator;
  
  edm::Handle<reco::CaloJetCollection> jets;
  
  MyRecHit hbheRecHit;
  MyRecHit hfRecHit;
  MyRecHit ebRecHit;
  MyRecHit eeRecHit;

  TTree* hbheTree;
  TTree* hfTree;
  TTree* ebTree;
  TTree* eeTree;
  double cone;

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

  edm::InputTag HcalRecHitHFSrc_;
  edm::InputTag HcalRecHitHBHESrc_;
  edm::InputTag EBSrc_;
  edm::InputTag EESrc_;
  edm::InputTag JetSrc_;
  bool excludeJets_;

};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
RecHitTreeProducer::RecHitTreeProducer(const edm::ParameterSet& iConfig) :
   cbins_(0),
   geo(0),
   cone(0.5)
{
   //now do what ever initialization is needed

  EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
  EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
  HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
  HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
  JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
  excludeJets_ = iConfig.getUntrackedParameter<bool>("excludeJets",false);
}


RecHitTreeProducer::~RecHitTreeProducer()
{
 
   // 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
RecHitTreeProducer::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
{

  hfRecHit.n = 0;
  hbheRecHit.n = 0;
  ebRecHit.n = 0;
  eeRecHit.n = 0;

  ev.getByLabel(EBSrc_,ebHits);
  ev.getByLabel(EESrc_,eeHits);

  ev.getByLabel(HcalRecHitHFSrc_,hfHits);
  ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);


  if(!excludeJets_) {
    ev.getByLabel(JetSrc_,jets);
  }
  
   if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);

   if(!geo){
      edm::ESHandle<CaloGeometry> pGeo;
      iSetup.get<CaloGeometryRecord>().get(pGeo);
      geo = pGeo.product();
   }
   
   for(unsigned int i = 0; i < hfHits->size(); ++i){
     const HFRecHit & hit= (*hfHits)[i];
     hfRecHit.e[hfRecHit.n] = hit.energy();
     hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
     hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
     hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
     hfRecHit.isjet[hfRecHit.n] = false;
     if(!excludeJets_){
       for(unsigned int j = 0 ; j < jets->size(); ++j){
         const reco::Jet& jet = (*jets)[j];
         double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
         if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
       }
     }
     hfRecHit.n++;
   }
   
   for(unsigned int i = 0; i < hbheHits->size(); ++i){
     const HBHERecHit & hit= (*hbheHits)[i];
     hbheRecHit.e[hbheRecHit.n] = hit.energy();
     hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
     hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
     hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
     hbheRecHit.isjet[hbheRecHit.n] = false; 
     if(!excludeJets_){
       for(unsigned int j = 0 ; j < jets->size(); ++j){
         const reco::Jet& jet = (*jets)[j];
         double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
         if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
       }
     }
     hbheRecHit.n++;
   }
   
   for(unsigned int i = 0; i < ebHits->size(); ++i){
     const EcalRecHit & hit= (*ebHits)[i];
     ebRecHit.e[ebRecHit.n] = hit.energy();
     ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
     ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
     ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
     ebRecHit.isjet[ebRecHit.n] = false;
     if(!excludeJets_){
       for(unsigned int j = 0 ; j < jets->size(); ++j){
         const reco::Jet& jet = (*jets)[j];
         double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
         if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
       }
     }
     ebRecHit.n++;
   }
   
   for(unsigned int i = 0; i < eeHits->size(); ++i){
     const EcalRecHit & hit= (*eeHits)[i];
     eeRecHit.e[eeRecHit.n] = hit.energy();
     eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
     eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
     eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
     eeRecHit.isjet[eeRecHit.n] = false;
     if(!excludeJets_){
       for(unsigned int j = 0 ; j < jets->size(); ++j){
         const reco::Jet& jet = (*jets)[j];
         double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
         if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
       }
     }
     eeRecHit.n++;
   }
   
   eeTree->Fill();
   ebTree->Fill();
   
   hbheTree->Fill();
   hfTree->Fill();
   
}


// ------------ method called once each job just before starting event loop  ------------
void 
RecHitTreeProducer::beginJob()
{
  
  hbheTree = fs->make<TTree>("hbhe","");
  hbheTree->Branch("n",&hbheRecHit.n,"n/I");
  hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
  hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
  hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
  hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
  hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/I");
  
  hfTree = fs->make<TTree>("hf","");
  hfTree->Branch("n",&hfRecHit.n,"n/I");
  hfTree->Branch("e",hfRecHit.e,"e[n]/F");
  hfTree->Branch("et",hfRecHit.et,"et[n]/F");
  hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
  hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
  hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/I");

  eeTree = fs->make<TTree>("ee","");
  eeTree->Branch("n",&eeRecHit.n,"n/I");
  eeTree->Branch("e",eeRecHit.e,"e[n]/F");
  eeTree->Branch("et",eeRecHit.et,"et[n]/F");
  eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
  eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
  eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/I");
 
  ebTree = fs->make<TTree>("eb","");
  ebTree->Branch("n",&ebRecHit.n,"n/I");
  ebTree->Branch("e",ebRecHit.e,"e[n]/F");
  ebTree->Branch("et",ebRecHit.et,"et[n]/F");
  ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
  ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
  ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");

}

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

double RecHitTreeProducer::getEt(const DetId &id, double energy){
  const GlobalPoint& pos=geo->getPosition(id);
  double et = energy*sin(pos.theta());
  return et;
}

double RecHitTreeProducer::getEta(const DetId &id){
  const GlobalPoint& pos=geo->getPosition(id);
  double et = pos.eta();
  return et;
}

double RecHitTreeProducer::getPhi(const DetId &id){
  const GlobalPoint& pos=geo->getPosition(id);
  double et = pos.phi();
  return et;
}


//define this as a plug-in
DEFINE_FWK_MODULE(RecHitTreeProducer);
Revision:	1.4
Committed:	Fri Oct 22 14:06:15 2010 UTC (14 years, 6 months ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Changes since 1.3:	+2 -2 lines
Log Message:	boolean branch
#	User	Rev	Content
1	yilmaz	1.1	// -- C++ --
2			//
3			// Package: RecHitTreeProducer
4			// Class: RecHitTreeProducer
5			//
6			/**\class RecHitTreeProducer RecHitTreeProducer.cc CmsHi/RecHitTreeProducer/src/RecHitTreeProducer.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.4	// $Id: RecHitTreeProducer.cc,v 1.3 2010/10/20 15:01:11 nart 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
56			#include "DataFormats/EcalRecHit/interface/EcalRecHitCollections.h"
57			#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
58			#include "DataFormats/HcalDetId/interface/HcalDetId.h"
59			#include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
60			#include "DataFormats/DetId/interface/DetId.h"
61
62
63			#include "TNtuple.h"
64
65			using namespace std;
66
67	nart	1.3	#define MAXHITS 1000000
68	yilmaz	1.1
69			struct MyRecHit{
70
71			int n;
72
73			float e[MAXHITS];
74			float et[MAXHITS];
75			float eta[MAXHITS];
76			float phi[MAXHITS];
77	nart	1.3	bool isjet[MAXHITS];
78	yilmaz	1.1
79			};
80
81
82
83			//
84			// class declaration
85			//
86
87			class RecHitTreeProducer : public edm::EDAnalyzer {
88			public:
89			explicit RecHitTreeProducer(const edm::ParameterSet&);
90			~RecHitTreeProducer();
91			double getEt(const DetId &id, double energy);
92			double getEta(const DetId &id);
93			double getPhi(const DetId &id);
94
95
96			private:
97			virtual void beginJob() ;
98			virtual void analyze(const edm::Event&, const edm::EventSetup&);
99			virtual void endJob() ;
100
101			// ----------member data ---------------------------
102
103			edm::Handle<reco::Centrality> cent;
104			edm::Handle<vector<double> > ktRhos;
105			edm::Handle<vector<double> > akRhos;
106
107	yilmaz	1.2	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits;
108			edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits;
109
110			edm::Handle<HFRecHitCollection> hfHits;
111			edm::Handle<HBHERecHitCollection> hbheHits;
112	yilmaz	1.1
113	nart	1.3	typedef vector<EcalRecHit>::const_iterator EcalIterator;
114
115			edm::Handle<reco::CaloJetCollection> jets;
116
117	yilmaz	1.1	MyRecHit hbheRecHit;
118			MyRecHit hfRecHit;
119	yilmaz	1.2	MyRecHit ebRecHit;
120			MyRecHit eeRecHit;
121	yilmaz	1.1
122			TTree* hbheTree;
123			TTree* hfTree;
124	yilmaz	1.2	TTree* ebTree;
125			TTree* eeTree;
126			double cone;
127	yilmaz	1.1
128			edm::Service<TFileService> fs;
129			const CentralityBins * cbins_;
130			const CaloGeometry *geo;
131
132			edm::InputTag HcalRecHitHFSrc_;
133			edm::InputTag HcalRecHitHBHESrc_;
134	yilmaz	1.2	edm::InputTag EBSrc_;
135			edm::InputTag EESrc_;
136	nart	1.3	edm::InputTag JetSrc_;
137			bool excludeJets_;
138	yilmaz	1.1
139			};
140
141			//
142			// constants, enums and typedefs
143			//
144
145			//
146			// static data member definitions
147			//
148
149			//
150			// constructors and destructor
151			//
152			RecHitTreeProducer::RecHitTreeProducer(const edm::ParameterSet& iConfig) :
153			cbins_(0),
154			geo(0),
155			cone(0.5)
156			{
157			//now do what ever initialization is needed
158	yilmaz	1.2
159			EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
160			EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
161	yilmaz	1.1	HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
162			HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
163	nart	1.3	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
164			excludeJets_ = iConfig.getUntrackedParameter<bool>("excludeJets",false);
165	yilmaz	1.1	}
166
167
168			RecHitTreeProducer::~RecHitTreeProducer()
169			{
170
171			// do anything here that needs to be done at desctruction time
172			// (e.g. close files, deallocate resources etc.)
173
174			}
175
176
177			//
178			// member functions
179			//
180
181			// ------------ method called to for each event ------------
182			void
183			RecHitTreeProducer::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
184			{
185
186			hfRecHit.n = 0;
187			hbheRecHit.n = 0;
188	yilmaz	1.2	ebRecHit.n = 0;
189			eeRecHit.n = 0;
190	yilmaz	1.1
191	yilmaz	1.2	ev.getByLabel(EBSrc_,ebHits);
192			ev.getByLabel(EESrc_,eeHits);
193	yilmaz	1.1
194			ev.getByLabel(HcalRecHitHFSrc_,hfHits);
195			ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
196
197	nart	1.3
198			if(!excludeJets_) {
199			ev.getByLabel(JetSrc_,jets);
200			}
201
202	yilmaz	1.1	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
203
204			if(!geo){
205			edm::ESHandle<CaloGeometry> pGeo;
206			iSetup.get<CaloGeometryRecord>().get(pGeo);
207			geo = pGeo.product();
208			}
209	nart	1.3
210	yilmaz	1.1	for(unsigned int i = 0; i < hfHits->size(); ++i){
211			const HFRecHit & hit= (*hfHits)[i];
212			hfRecHit.e[hfRecHit.n] = hit.energy();
213			hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
214			hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
215			hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
216	nart	1.3	hfRecHit.isjet[hfRecHit.n] = false;
217			if(!excludeJets_){
218			for(unsigned int j = 0 ; j < jets->size(); ++j){
219			const reco::Jet& jet = (*jets)[j];
220			double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
221			if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
222			}
223			}
224	yilmaz	1.1	hfRecHit.n++;
225			}
226	nart	1.3
227	yilmaz	1.1	for(unsigned int i = 0; i < hbheHits->size(); ++i){
228			const HBHERecHit & hit= (*hbheHits)[i];
229	nart	1.3	hbheRecHit.e[hbheRecHit.n] = hit.energy();
230			hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
231			hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
232			hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
233			hbheRecHit.isjet[hbheRecHit.n] = false;
234			if(!excludeJets_){
235			for(unsigned int j = 0 ; j < jets->size(); ++j){
236			const reco::Jet& jet = (*jets)[j];
237			double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
238			if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
239			}
240			}
241	yilmaz	1.2	hbheRecHit.n++;
242			}
243	nart	1.3
244	yilmaz	1.2	for(unsigned int i = 0; i < ebHits->size(); ++i){
245			const EcalRecHit & hit= (*ebHits)[i];
246			ebRecHit.e[ebRecHit.n] = hit.energy();
247			ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
248			ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
249			ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
250	nart	1.3	ebRecHit.isjet[ebRecHit.n] = false;
251			if(!excludeJets_){
252			for(unsigned int j = 0 ; j < jets->size(); ++j){
253			const reco::Jet& jet = (*jets)[j];
254			double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
255			if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
256			}
257			}
258	yilmaz	1.2	ebRecHit.n++;
259			}
260	nart	1.3
261	yilmaz	1.2	for(unsigned int i = 0; i < eeHits->size(); ++i){
262			const EcalRecHit & hit= (*eeHits)[i];
263			eeRecHit.e[eeRecHit.n] = hit.energy();
264			eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
265			eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
266			eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
267	nart	1.3	eeRecHit.isjet[eeRecHit.n] = false;
268			if(!excludeJets_){
269			for(unsigned int j = 0 ; j < jets->size(); ++j){
270			const reco::Jet& jet = (*jets)[j];
271			double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
272			if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
273			}
274			}
275	yilmaz	1.2	eeRecHit.n++;
276	yilmaz	1.1	}
277	nart	1.3
278	yilmaz	1.2	eeTree->Fill();
279			ebTree->Fill();
280	nart	1.3
281	yilmaz	1.1	hbheTree->Fill();
282			hfTree->Fill();
283	nart	1.3
284	yilmaz	1.1	}
285
286
287			// ------------ method called once each job just before starting event loop ------------
288			void
289			RecHitTreeProducer::beginJob()
290			{
291	nart	1.3
292	yilmaz	1.1	hbheTree = fs->make<TTree>("hbhe","");
293			hbheTree->Branch("n",&hbheRecHit.n,"n/I");
294			hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
295			hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
296			hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
297			hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
298	nart	1.3	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/I");
299
300	yilmaz	1.1	hfTree = fs->make<TTree>("hf","");
301			hfTree->Branch("n",&hfRecHit.n,"n/I");
302			hfTree->Branch("e",hfRecHit.e,"e[n]/F");
303			hfTree->Branch("et",hfRecHit.et,"et[n]/F");
304			hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
305			hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
306	nart	1.3	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/I");
307	yilmaz	1.1
308	yilmaz	1.2	eeTree = fs->make<TTree>("ee","");
309			eeTree->Branch("n",&eeRecHit.n,"n/I");
310			eeTree->Branch("e",eeRecHit.e,"e[n]/F");
311			eeTree->Branch("et",eeRecHit.et,"et[n]/F");
312			eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
313			eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
314	nart	1.3	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/I");
315
316	yilmaz	1.2	ebTree = fs->make<TTree>("eb","");
317			ebTree->Branch("n",&ebRecHit.n,"n/I");
318			ebTree->Branch("e",ebRecHit.e,"e[n]/F");
319			ebTree->Branch("et",ebRecHit.et,"et[n]/F");
320			ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
321			ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
322	yilmaz	1.4	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
323	yilmaz	1.2
324	yilmaz	1.1	}
325
326			// ------------ method called once each job just after ending the event loop ------------
327			void
328			RecHitTreeProducer::endJob() {
329			}
330
331			double RecHitTreeProducer::getEt(const DetId &id, double energy){
332			const GlobalPoint& pos=geo->getPosition(id);
333			double et = energy*sin(pos.theta());
334			return et;
335			}
336	nart	1.3
337	yilmaz	1.1	double RecHitTreeProducer::getEta(const DetId &id){
338			const GlobalPoint& pos=geo->getPosition(id);
339			double et = pos.eta();
340			return et;
341			}
342
343			double RecHitTreeProducer::getPhi(const DetId &id){
344			const GlobalPoint& pos=geo->getPosition(id);
345			double et = pos.phi();
346			return et;
347			}
348
349
350
351			//define this as a plug-in
352			DEFINE_FWK_MODULE(RecHitTreeProducer);