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.5 2010/10/22 14:08:05 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/EcalRecHit/interface/EcalRecHitCollections.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/HcalDetId/interface/HcalDetId.h"
#include "DataFormats/EcalDetId/interface/EcalDetIdCollections.h"
#include "DataFormats/EgammaReco/interface/BasicClusterFwd.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];

   float jtpt;
   float jteta;
   float jtphi;

};


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

   edm::Handle<reco::BasicClusterCollection> bClusters;
   edm::Handle<CaloTowerCollection> towers;


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

  TTree* hbheTree;
  TTree* hfTree;
  TTree* ebTree;
  TTree* eeTree;
   TTree* bcTree;
   TTree* towerTree;

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

   edm::InputTag TowerSrc_;

  edm::InputTag JetSrc_;

  bool useJets_;
   bool doBasicClusters_;

};

//
// 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"));
  BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
  TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
  JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
  useJets_ = iConfig.getUntrackedParameter<bool>("useJets",false);
  doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",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(useJets_) {
    ev.getByLabel(JetSrc_,jets);
  }

  if(doBasicClusters_){
     ev.getByLabel(BCSrc_,bClusters);
  }


  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(!useJets_){
       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(!useJets_){
       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(!useJets_){
       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(!useJets_){
       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++;
   }


   if(doBasicClusters_){
      for(unsigned int j = 0 ; j < jets->size(); ++j){
         const reco::Jet& jet = (*jets)[j];
         myBC.n = 0;
         myBC.jtpt = jet.pt();
         myBC.jteta = jet.eta();
         myBC.jtphi = jet.phi();

         for(unsigned int i = 0; i < bClusters->size(); ++i){
            const reco::BasicCluster & bc= (*bClusters)[i];
            double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
            if(dr < cone){ 
               myBC.e[myBC.n] = bc.energy();
               myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
               myBC.eta[myBC.n] = bc.eta();
               myBC.phi[myBC.n] = bc.phi();
               myBC.n++; 
            }
         }
         bcTree->Fill(); 
      }
   }
  
   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]/O");
  
  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]/O");

  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]/O");
 
  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");

  if(doBasicClusters_){
     bcTree = fs->make<TTree>("bc","");
     bcTree->Branch("n",&myBC.n,"n/I");
     bcTree->Branch("e",myBC.e,"e[n]/F");
     bcTree->Branch("et",myBC.et,"et[n]/F");
     bcTree->Branch("eta",myBC.eta,"eta[n]/F");
     bcTree->Branch("phi",myBC.phi,"phi[n]/F");
     bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
     bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
     bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
     //     bcTree->Branch("isjet",bcRecHit.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.6
Committed:	Mon Nov 1 21:48:31 2010 UTC (14 years, 6 months ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Changes since 1.5:	+79 -12 lines
Log Message:	more info: basic clusters, towers
#	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.6	// $Id: RecHitTreeProducer.cc,v 1.5 2010/10/22 14:08:05 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
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	yilmaz	1.6	#include "DataFormats/EgammaReco/interface/BasicClusterFwd.h"
61	yilmaz	1.1	#include "DataFormats/DetId/interface/DetId.h"
62
63
64			#include "TNtuple.h"
65
66			using namespace std;
67
68	nart	1.3	#define MAXHITS 1000000
69	yilmaz	1.1
70			struct MyRecHit{
71
72			int n;
73
74			float e[MAXHITS];
75			float et[MAXHITS];
76			float eta[MAXHITS];
77			float phi[MAXHITS];
78	nart	1.3	bool isjet[MAXHITS];
79	yilmaz	1.1
80	yilmaz	1.6	float jtpt;
81			float jteta;
82			float jtphi;
83
84	yilmaz	1.1	};
85
86
87
88			//
89			// class declaration
90			//
91
92			class RecHitTreeProducer : public edm::EDAnalyzer {
93			public:
94			explicit RecHitTreeProducer(const edm::ParameterSet&);
95			~RecHitTreeProducer();
96			double getEt(const DetId &id, double energy);
97			double getEta(const DetId &id);
98			double getPhi(const DetId &id);
99
100
101			private:
102			virtual void beginJob() ;
103			virtual void analyze(const edm::Event&, const edm::EventSetup&);
104			virtual void endJob() ;
105
106			// ----------member data ---------------------------
107
108			edm::Handle<reco::Centrality> cent;
109			edm::Handle<vector<double> > ktRhos;
110			edm::Handle<vector<double> > akRhos;
111
112	yilmaz	1.2	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits;
113			edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits;
114
115			edm::Handle<HFRecHitCollection> hfHits;
116			edm::Handle<HBHERecHitCollection> hbheHits;
117	yilmaz	1.1
118	yilmaz	1.6	edm::Handle<reco::BasicClusterCollection> bClusters;
119			edm::Handle<CaloTowerCollection> towers;
120
121
122	nart	1.3	typedef vector<EcalRecHit>::const_iterator EcalIterator;
123
124			edm::Handle<reco::CaloJetCollection> jets;
125
126	yilmaz	1.1	MyRecHit hbheRecHit;
127			MyRecHit hfRecHit;
128	yilmaz	1.2	MyRecHit ebRecHit;
129			MyRecHit eeRecHit;
130	yilmaz	1.6	MyRecHit myBC;
131			MyRecHit myTowers;
132	yilmaz	1.1
133			TTree* hbheTree;
134			TTree* hfTree;
135	yilmaz	1.2	TTree* ebTree;
136			TTree* eeTree;
137	yilmaz	1.6	TTree* bcTree;
138			TTree* towerTree;
139
140	yilmaz	1.2	double cone;
141	yilmaz	1.1
142			edm::Service<TFileService> fs;
143			const CentralityBins * cbins_;
144			const CaloGeometry *geo;
145
146			edm::InputTag HcalRecHitHFSrc_;
147			edm::InputTag HcalRecHitHBHESrc_;
148	yilmaz	1.2	edm::InputTag EBSrc_;
149			edm::InputTag EESrc_;
150	yilmaz	1.6	edm::InputTag BCSrc_;
151
152			edm::InputTag TowerSrc_;
153
154	nart	1.3	edm::InputTag JetSrc_;
155	yilmaz	1.6
156			bool useJets_;
157			bool doBasicClusters_;
158	yilmaz	1.1
159			};
160
161			//
162			// constants, enums and typedefs
163			//
164
165			//
166			// static data member definitions
167			//
168
169			//
170			// constructors and destructor
171			//
172			RecHitTreeProducer::RecHitTreeProducer(const edm::ParameterSet& iConfig) :
173			cbins_(0),
174			geo(0),
175			cone(0.5)
176			{
177			//now do what ever initialization is needed
178	yilmaz	1.2	EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
179			EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
180	yilmaz	1.1	HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
181			HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
182	yilmaz	1.6	BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
183			TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
184	nart	1.3	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
185	yilmaz	1.6	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",false);
186			doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
187
188	yilmaz	1.1	}
189
190
191			RecHitTreeProducer::~RecHitTreeProducer()
192			{
193
194			// do anything here that needs to be done at desctruction time
195			// (e.g. close files, deallocate resources etc.)
196
197			}
198
199
200			//
201			// member functions
202			//
203
204			// ------------ method called to for each event ------------
205			void
206			RecHitTreeProducer::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
207			{
208
209			hfRecHit.n = 0;
210			hbheRecHit.n = 0;
211	yilmaz	1.2	ebRecHit.n = 0;
212			eeRecHit.n = 0;
213	yilmaz	1.1
214	yilmaz	1.2	ev.getByLabel(EBSrc_,ebHits);
215			ev.getByLabel(EESrc_,eeHits);
216	yilmaz	1.1
217			ev.getByLabel(HcalRecHitHFSrc_,hfHits);
218			ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
219
220	yilmaz	1.6	if(useJets_) {
221			ev.getByLabel(JetSrc_,jets);
222			}
223	nart	1.3
224	yilmaz	1.6	if(doBasicClusters_){
225			ev.getByLabel(BCSrc_,bClusters);
226	nart	1.3	}
227	yilmaz	1.6
228
229	nart	1.3
230	yilmaz	1.6	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
231	yilmaz	1.1
232			if(!geo){
233			edm::ESHandle<CaloGeometry> pGeo;
234			iSetup.get<CaloGeometryRecord>().get(pGeo);
235			geo = pGeo.product();
236			}
237	nart	1.3
238	yilmaz	1.1	for(unsigned int i = 0; i < hfHits->size(); ++i){
239			const HFRecHit & hit= (*hfHits)[i];
240			hfRecHit.e[hfRecHit.n] = hit.energy();
241			hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
242			hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
243			hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
244	nart	1.3	hfRecHit.isjet[hfRecHit.n] = false;
245	yilmaz	1.6	if(!useJets_){
246	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
247			const reco::Jet& jet = (*jets)[j];
248			double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
249			if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
250			}
251			}
252	yilmaz	1.1	hfRecHit.n++;
253			}
254	nart	1.3
255	yilmaz	1.1	for(unsigned int i = 0; i < hbheHits->size(); ++i){
256			const HBHERecHit & hit= (*hbheHits)[i];
257	nart	1.3	hbheRecHit.e[hbheRecHit.n] = hit.energy();
258			hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
259			hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
260			hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
261			hbheRecHit.isjet[hbheRecHit.n] = false;
262	yilmaz	1.6	if(!useJets_){
263	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
264			const reco::Jet& jet = (*jets)[j];
265			double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
266			if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
267			}
268			}
269	yilmaz	1.2	hbheRecHit.n++;
270			}
271	nart	1.3
272	yilmaz	1.2	for(unsigned int i = 0; i < ebHits->size(); ++i){
273			const EcalRecHit & hit= (*ebHits)[i];
274			ebRecHit.e[ebRecHit.n] = hit.energy();
275			ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
276			ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
277			ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
278	nart	1.3	ebRecHit.isjet[ebRecHit.n] = false;
279	yilmaz	1.6	if(!useJets_){
280	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
281			const reco::Jet& jet = (*jets)[j];
282			double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
283			if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
284			}
285			}
286	yilmaz	1.2	ebRecHit.n++;
287			}
288	nart	1.3
289	yilmaz	1.2	for(unsigned int i = 0; i < eeHits->size(); ++i){
290			const EcalRecHit & hit= (*eeHits)[i];
291			eeRecHit.e[eeRecHit.n] = hit.energy();
292			eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
293			eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
294			eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
295	nart	1.3	eeRecHit.isjet[eeRecHit.n] = false;
296	yilmaz	1.6	if(!useJets_){
297	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
298			const reco::Jet& jet = (*jets)[j];
299			double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
300			if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
301			}
302			}
303	yilmaz	1.2	eeRecHit.n++;
304	yilmaz	1.1	}
305	yilmaz	1.6
306
307			if(doBasicClusters_){
308			for(unsigned int j = 0 ; j < jets->size(); ++j){
309			const reco::Jet& jet = (*jets)[j];
310			myBC.n = 0;
311			myBC.jtpt = jet.pt();
312			myBC.jteta = jet.eta();
313			myBC.jtphi = jet.phi();
314
315			for(unsigned int i = 0; i < bClusters->size(); ++i){
316			const reco::BasicCluster & bc= (*bClusters)[i];
317			double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
318			if(dr < cone){
319			myBC.e[myBC.n] = bc.energy();
320			myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
321			myBC.eta[myBC.n] = bc.eta();
322			myBC.phi[myBC.n] = bc.phi();
323			myBC.n++;
324			}
325			}
326			bcTree->Fill();
327			}
328			}
329
330	yilmaz	1.2	eeTree->Fill();
331			ebTree->Fill();
332	nart	1.3
333	yilmaz	1.1	hbheTree->Fill();
334			hfTree->Fill();
335	nart	1.3
336	yilmaz	1.1	}
337
338
339			// ------------ method called once each job just before starting event loop ------------
340			void
341			RecHitTreeProducer::beginJob()
342			{
343	nart	1.3
344	yilmaz	1.1	hbheTree = fs->make<TTree>("hbhe","");
345			hbheTree->Branch("n",&hbheRecHit.n,"n/I");
346			hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
347			hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
348			hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
349			hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
350	yilmaz	1.5	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
351	nart	1.3
352	yilmaz	1.1	hfTree = fs->make<TTree>("hf","");
353			hfTree->Branch("n",&hfRecHit.n,"n/I");
354			hfTree->Branch("e",hfRecHit.e,"e[n]/F");
355			hfTree->Branch("et",hfRecHit.et,"et[n]/F");
356			hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
357			hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
358	yilmaz	1.5	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
359	yilmaz	1.1
360	yilmaz	1.2	eeTree = fs->make<TTree>("ee","");
361			eeTree->Branch("n",&eeRecHit.n,"n/I");
362			eeTree->Branch("e",eeRecHit.e,"e[n]/F");
363			eeTree->Branch("et",eeRecHit.et,"et[n]/F");
364			eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
365			eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
366	yilmaz	1.5	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
367	nart	1.3
368	yilmaz	1.2	ebTree = fs->make<TTree>("eb","");
369			ebTree->Branch("n",&ebRecHit.n,"n/I");
370			ebTree->Branch("e",ebRecHit.e,"e[n]/F");
371			ebTree->Branch("et",ebRecHit.et,"et[n]/F");
372			ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
373			ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
374	yilmaz	1.4	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
375	yilmaz	1.2
376	yilmaz	1.6	if(doBasicClusters_){
377			bcTree = fs->make<TTree>("bc","");
378			bcTree->Branch("n",&myBC.n,"n/I");
379			bcTree->Branch("e",myBC.e,"e[n]/F");
380			bcTree->Branch("et",myBC.et,"et[n]/F");
381			bcTree->Branch("eta",myBC.eta,"eta[n]/F");
382			bcTree->Branch("phi",myBC.phi,"phi[n]/F");
383			bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
384			bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
385			bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
386			// bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
387			}
388
389
390
391	yilmaz	1.1	}
392
393			// ------------ method called once each job just after ending the event loop ------------
394			void
395			RecHitTreeProducer::endJob() {
396			}
397
398			double RecHitTreeProducer::getEt(const DetId &id, double energy){
399			const GlobalPoint& pos=geo->getPosition(id);
400			double et = energy*sin(pos.theta());
401			return et;
402			}
403	nart	1.3
404	yilmaz	1.1	double RecHitTreeProducer::getEta(const DetId &id){
405			const GlobalPoint& pos=geo->getPosition(id);
406			double et = pos.eta();
407			return et;
408			}
409
410			double RecHitTreeProducer::getPhi(const DetId &id){
411			const GlobalPoint& pos=geo->getPosition(id);
412			double et = pos.phi();
413			return et;
414			}
415
416
417
418			//define this as a plug-in
419			DEFINE_FWK_MODULE(RecHitTreeProducer);