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.6 2010/11/01 21:48:31 yilmaz Exp $
//
//


// system include files
#include <memory>
#include <vector>
#include <iostream>

// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/DetId/interface/DetId.h"
#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
#include "Geometry/Records/interface/IdealGeometryRecord.h"
#include "Geometry/Records/interface/CaloGeometryRecord.h"
#include "DataFormats/GeometryVector/interface/GlobalPoint.h"

#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "FWCore/ServiceRegistry/interface/Service.h"

#include "DataFormats/Math/interface/deltaR.h"

#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/FWLite/interface/ChainEvent.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
#include "DataFormats/CaloTowers/interface/CaloTower.h"
#include "DataFormats/HeavyIonEvent/interface/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;

};


struct MyBkg{
   int n;
   float rho[50];
   float sigma[50];
};


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

   edm::Handle<std::vector<double> > rhos;
   edm::Handle<std::vector<double> > sigmas;
  
  MyRecHit hbheRecHit;
  MyRecHit hfRecHit;
  MyRecHit ebRecHit;
  MyRecHit eeRecHit;
   MyRecHit myBC;
   MyRecHit myTowers;
   MyBkg bkg;


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

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

   edm::InputTag FastJetTag_;

  bool useJets_;
   bool doBasicClusters_;
   bool doTowers_;
   bool doEcal_;
   bool doHcal_;

   bool doFastJet_;
};

//
// 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",true);
  doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
  doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
  doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
  doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
  doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
  FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));

}


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;
  myBC.n = 0;
   myTowers.n = 0;
   bkg.n = 0;

   if(doEcal_){
  ev.getByLabel(EBSrc_,ebHits);
  ev.getByLabel(EESrc_,eeHits);
   }
  if(doHcal_){
  ev.getByLabel(HcalRecHitHFSrc_,hfHits);
  ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
  }
  if(useJets_) {
    ev.getByLabel(JetSrc_,jets);
  }

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

  if(doTowers_){
     ev.getByLabel(TowerSrc_,towers);
  }

  if(doFastJet_){
     ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
     ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
     bkg.n = rhos->size();
     for(unsigned int i = 0; i < rhos->size(); ++i){
        bkg.rho[i] = (*rhos)[i];
        bkg.sigma[i] = (*sigmas)[i];
     }
  }
  
  if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);

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

      for(unsigned int i = 0; i < towers->size(); ++i){
      const CaloTower & hit= (*towers)[i];
      myTowers.e[myTowers.n] = hit.energy();
      myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
      myTowers.eta[myTowers.n] = getEta(hit.id());
      myTowers.phi[myTowers.n] = getPhi(hit.id());
      myTowers.isjet[myTowers.n] = false;
      if(useJets_){
         for(unsigned int j = 0 ; j < jets->size(); ++j){
            const reco::Jet& jet = (*jets)[j];
            double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
            if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
         }
      }
      myTowers.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(); 
      }
   }

   towerTree->Fill();
  
   eeTree->Fill();
   ebTree->Fill();
   
   hbheTree->Fill();
   hfTree->Fill();
   bkgTree->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");

  towerTree = fs->make<TTree>("tower","");
  towerTree->Branch("n",&myTowers.n,"n/I");
  towerTree->Branch("e",myTowers.e,"e[n]/F");
  towerTree->Branch("et",myTowers.et,"et[n]/F");
  towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
  towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
  towerTree->Branch("isjet",myTowers.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");
  }

  if(doFastJet_){
     bkgTree = fs->make<TTree>("bkg","");
     bkgTree->Branch("n",&bkg.n,"n/I");
     bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
     bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
  }

}

// ------------ 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.7
Committed:	Tue Nov 2 12:58:33 2010 UTC (14 years, 6 months ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Changes since 1.6:	+90 -12 lines
Log Message:	update
#	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.7	// $Id: RecHitTreeProducer.cc,v 1.6 2010/11/01 21:48:31 yilmaz Exp $
17	yilmaz	1.1	//
18			//
19
20
21			// system include files
22			#include <memory>
23			#include <vector>
24			#include <iostream>
25
26			// user include files
27			#include "FWCore/Framework/interface/Frameworkfwd.h"
28			#include "FWCore/Framework/interface/EDAnalyzer.h"
29
30			#include "FWCore/Framework/interface/Event.h"
31			#include "FWCore/Framework/interface/MakerMacros.h"
32			#include "FWCore/Framework/interface/ESHandle.h"
33			#include "FWCore/ParameterSet/interface/ParameterSet.h"
34
35			#include "DataFormats/DetId/interface/DetId.h"
36			#include "Geometry/CaloGeometry/interface/CaloGeometry.h"
37			#include "Geometry/Records/interface/IdealGeometryRecord.h"
38			#include "Geometry/Records/interface/CaloGeometryRecord.h"
39			#include "DataFormats/GeometryVector/interface/GlobalPoint.h"
40
41			#include "CommonTools/UtilAlgos/interface/TFileService.h"
42			#include "FWCore/ServiceRegistry/interface/Service.h"
43
44			#include "DataFormats/Math/interface/deltaR.h"
45
46			#include "DataFormats/Common/interface/Handle.h"
47			#include "DataFormats/FWLite/interface/ChainEvent.h"
48			#include "FWCore/Utilities/interface/InputTag.h"
49			#include "DataFormats/HeavyIonEvent/interface/CentralityBins.h"
50			#include "DataFormats/CaloTowers/interface/CaloTower.h"
51			#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	yilmaz	1.7	struct MyBkg{
88			int n;
89			float rho[50];
90			float sigma[50];
91			};
92
93	yilmaz	1.1
94			//
95			// class declaration
96			//
97
98			class RecHitTreeProducer : public edm::EDAnalyzer {
99			public:
100			explicit RecHitTreeProducer(const edm::ParameterSet&);
101			~RecHitTreeProducer();
102			double getEt(const DetId &id, double energy);
103			double getEta(const DetId &id);
104			double getPhi(const DetId &id);
105
106
107			private:
108			virtual void beginJob() ;
109			virtual void analyze(const edm::Event&, const edm::EventSetup&);
110			virtual void endJob() ;
111
112			// ----------member data ---------------------------
113
114			edm::Handle<reco::Centrality> cent;
115			edm::Handle<vector<double> > ktRhos;
116			edm::Handle<vector<double> > akRhos;
117
118	yilmaz	1.2	edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > ebHits;
119			edm::Handle<edm::SortedCollection<EcalRecHit,edm::StrictWeakOrdering<EcalRecHit> > > eeHits;
120
121			edm::Handle<HFRecHitCollection> hfHits;
122			edm::Handle<HBHERecHitCollection> hbheHits;
123	yilmaz	1.1
124	yilmaz	1.6	edm::Handle<reco::BasicClusterCollection> bClusters;
125			edm::Handle<CaloTowerCollection> towers;
126
127
128	nart	1.3	typedef vector<EcalRecHit>::const_iterator EcalIterator;
129
130			edm::Handle<reco::CaloJetCollection> jets;
131	yilmaz	1.7
132			edm::Handle<std::vector<double> > rhos;
133			edm::Handle<std::vector<double> > sigmas;
134	nart	1.3
135	yilmaz	1.1	MyRecHit hbheRecHit;
136			MyRecHit hfRecHit;
137	yilmaz	1.2	MyRecHit ebRecHit;
138			MyRecHit eeRecHit;
139	yilmaz	1.6	MyRecHit myBC;
140			MyRecHit myTowers;
141	yilmaz	1.7	MyBkg bkg;
142
143	yilmaz	1.1
144			TTree* hbheTree;
145			TTree* hfTree;
146	yilmaz	1.2	TTree* ebTree;
147			TTree* eeTree;
148	yilmaz	1.6	TTree* bcTree;
149			TTree* towerTree;
150	yilmaz	1.7	TTree* bkgTree;
151	yilmaz	1.6
152	yilmaz	1.2	double cone;
153	yilmaz	1.1
154			edm::Service<TFileService> fs;
155			const CentralityBins * cbins_;
156			const CaloGeometry *geo;
157
158			edm::InputTag HcalRecHitHFSrc_;
159			edm::InputTag HcalRecHitHBHESrc_;
160	yilmaz	1.2	edm::InputTag EBSrc_;
161			edm::InputTag EESrc_;
162	yilmaz	1.6	edm::InputTag BCSrc_;
163			edm::InputTag TowerSrc_;
164
165	nart	1.3	edm::InputTag JetSrc_;
166	yilmaz	1.6
167	yilmaz	1.7	edm::InputTag FastJetTag_;
168
169	yilmaz	1.6	bool useJets_;
170			bool doBasicClusters_;
171	yilmaz	1.7	bool doTowers_;
172			bool doEcal_;
173			bool doHcal_;
174	yilmaz	1.1
175	yilmaz	1.7	bool doFastJet_;
176	yilmaz	1.1	};
177
178			//
179			// constants, enums and typedefs
180			//
181
182			//
183			// static data member definitions
184			//
185
186			//
187			// constructors and destructor
188			//
189			RecHitTreeProducer::RecHitTreeProducer(const edm::ParameterSet& iConfig) :
190			cbins_(0),
191			geo(0),
192			cone(0.5)
193			{
194			//now do what ever initialization is needed
195	yilmaz	1.2	EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
196			EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
197	yilmaz	1.1	HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
198			HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
199	yilmaz	1.6	BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
200			TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
201	nart	1.3	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeCone5CaloJets"));
202	yilmaz	1.7	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
203	yilmaz	1.6	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
204	yilmaz	1.7	doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
205			doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
206			doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
207			doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
208			FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
209	yilmaz	1.6
210	yilmaz	1.1	}
211
212
213			RecHitTreeProducer::~RecHitTreeProducer()
214			{
215
216			// do anything here that needs to be done at desctruction time
217			// (e.g. close files, deallocate resources etc.)
218
219			}
220
221
222			//
223			// member functions
224			//
225
226			// ------------ method called to for each event ------------
227			void
228			RecHitTreeProducer::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
229			{
230
231			hfRecHit.n = 0;
232			hbheRecHit.n = 0;
233	yilmaz	1.2	ebRecHit.n = 0;
234			eeRecHit.n = 0;
235	yilmaz	1.7	myBC.n = 0;
236			myTowers.n = 0;
237			bkg.n = 0;
238	yilmaz	1.1
239	yilmaz	1.7	if(doEcal_){
240	yilmaz	1.2	ev.getByLabel(EBSrc_,ebHits);
241			ev.getByLabel(EESrc_,eeHits);
242	yilmaz	1.7	}
243			if(doHcal_){
244	yilmaz	1.1	ev.getByLabel(HcalRecHitHFSrc_,hfHits);
245			ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
246	yilmaz	1.7	}
247	yilmaz	1.6	if(useJets_) {
248			ev.getByLabel(JetSrc_,jets);
249			}
250	nart	1.3
251	yilmaz	1.6	if(doBasicClusters_){
252			ev.getByLabel(BCSrc_,bClusters);
253	nart	1.3	}
254	yilmaz	1.6
255	yilmaz	1.7	if(doTowers_){
256			ev.getByLabel(TowerSrc_,towers);
257			}
258	yilmaz	1.6
259	yilmaz	1.7	if(doFastJet_){
260			ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
261			ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
262			bkg.n = rhos->size();
263			for(unsigned int i = 0; i < rhos->size(); ++i){
264			bkg.rho[i] = (*rhos)[i];
265			bkg.sigma[i] = (*sigmas)[i];
266			}
267			}
268	nart	1.3
269	yilmaz	1.6	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
270	yilmaz	1.1
271			if(!geo){
272			edm::ESHandle<CaloGeometry> pGeo;
273			iSetup.get<CaloGeometryRecord>().get(pGeo);
274			geo = pGeo.product();
275			}
276	nart	1.3
277	yilmaz	1.7	if(doHcal_){
278	yilmaz	1.1	for(unsigned int i = 0; i < hfHits->size(); ++i){
279			const HFRecHit & hit= (*hfHits)[i];
280			hfRecHit.e[hfRecHit.n] = hit.energy();
281			hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
282			hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
283			hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
284	nart	1.3	hfRecHit.isjet[hfRecHit.n] = false;
285	yilmaz	1.7	if(useJets_){
286	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
287			const reco::Jet& jet = (*jets)[j];
288			double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
289			if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
290			}
291			}
292	yilmaz	1.1	hfRecHit.n++;
293			}
294	nart	1.3
295	yilmaz	1.1	for(unsigned int i = 0; i < hbheHits->size(); ++i){
296			const HBHERecHit & hit= (*hbheHits)[i];
297	nart	1.3	hbheRecHit.e[hbheRecHit.n] = hit.energy();
298			hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
299			hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
300			hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
301			hbheRecHit.isjet[hbheRecHit.n] = false;
302	yilmaz	1.7	if(useJets_){
303	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
304			const reco::Jet& jet = (*jets)[j];
305			double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
306			if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
307			}
308			}
309	yilmaz	1.2	hbheRecHit.n++;
310			}
311	yilmaz	1.7	}
312			if(doEcal_){
313	yilmaz	1.2	for(unsigned int i = 0; i < ebHits->size(); ++i){
314			const EcalRecHit & hit= (*ebHits)[i];
315			ebRecHit.e[ebRecHit.n] = hit.energy();
316			ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
317			ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
318			ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
319	nart	1.3	ebRecHit.isjet[ebRecHit.n] = false;
320	yilmaz	1.7	if(useJets_){
321	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
322			const reco::Jet& jet = (*jets)[j];
323			double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
324			if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
325			}
326			}
327	yilmaz	1.2	ebRecHit.n++;
328			}
329	nart	1.3
330	yilmaz	1.2	for(unsigned int i = 0; i < eeHits->size(); ++i){
331			const EcalRecHit & hit= (*eeHits)[i];
332			eeRecHit.e[eeRecHit.n] = hit.energy();
333			eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
334			eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
335			eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
336	nart	1.3	eeRecHit.isjet[eeRecHit.n] = false;
337	yilmaz	1.7	if(useJets_){
338	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
339			const reco::Jet& jet = (*jets)[j];
340			double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
341			if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
342			}
343			}
344	yilmaz	1.2	eeRecHit.n++;
345	yilmaz	1.1	}
346	yilmaz	1.7	}
347
348			if(doTowers_){
349
350			for(unsigned int i = 0; i < towers->size(); ++i){
351			const CaloTower & hit= (*towers)[i];
352			myTowers.e[myTowers.n] = hit.energy();
353			myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
354			myTowers.eta[myTowers.n] = getEta(hit.id());
355			myTowers.phi[myTowers.n] = getPhi(hit.id());
356			myTowers.isjet[myTowers.n] = false;
357			if(useJets_){
358			for(unsigned int j = 0 ; j < jets->size(); ++j){
359			const reco::Jet& jet = (*jets)[j];
360			double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
361			if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
362			}
363			}
364			myTowers.n++;
365			}
366	yilmaz	1.6
367	yilmaz	1.7	}
368	yilmaz	1.6
369			if(doBasicClusters_){
370			for(unsigned int j = 0 ; j < jets->size(); ++j){
371			const reco::Jet& jet = (*jets)[j];
372			myBC.n = 0;
373			myBC.jtpt = jet.pt();
374			myBC.jteta = jet.eta();
375			myBC.jtphi = jet.phi();
376
377			for(unsigned int i = 0; i < bClusters->size(); ++i){
378			const reco::BasicCluster & bc= (*bClusters)[i];
379			double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
380			if(dr < cone){
381			myBC.e[myBC.n] = bc.energy();
382			myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
383			myBC.eta[myBC.n] = bc.eta();
384			myBC.phi[myBC.n] = bc.phi();
385			myBC.n++;
386			}
387			}
388			bcTree->Fill();
389			}
390			}
391	yilmaz	1.7
392			towerTree->Fill();
393	yilmaz	1.6
394	yilmaz	1.2	eeTree->Fill();
395			ebTree->Fill();
396	nart	1.3
397	yilmaz	1.1	hbheTree->Fill();
398			hfTree->Fill();
399	yilmaz	1.7	bkgTree->Fill();
400	yilmaz	1.1	}
401
402
403			// ------------ method called once each job just before starting event loop ------------
404			void
405			RecHitTreeProducer::beginJob()
406			{
407	nart	1.3
408	yilmaz	1.1	hbheTree = fs->make<TTree>("hbhe","");
409			hbheTree->Branch("n",&hbheRecHit.n,"n/I");
410			hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
411			hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
412			hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
413			hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
414	yilmaz	1.5	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
415	nart	1.3
416	yilmaz	1.1	hfTree = fs->make<TTree>("hf","");
417			hfTree->Branch("n",&hfRecHit.n,"n/I");
418			hfTree->Branch("e",hfRecHit.e,"e[n]/F");
419			hfTree->Branch("et",hfRecHit.et,"et[n]/F");
420			hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
421			hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
422	yilmaz	1.5	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
423	yilmaz	1.1
424	yilmaz	1.2	eeTree = fs->make<TTree>("ee","");
425			eeTree->Branch("n",&eeRecHit.n,"n/I");
426			eeTree->Branch("e",eeRecHit.e,"e[n]/F");
427			eeTree->Branch("et",eeRecHit.et,"et[n]/F");
428			eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
429			eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
430	yilmaz	1.5	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
431	nart	1.3
432	yilmaz	1.2	ebTree = fs->make<TTree>("eb","");
433			ebTree->Branch("n",&ebRecHit.n,"n/I");
434			ebTree->Branch("e",ebRecHit.e,"e[n]/F");
435			ebTree->Branch("et",ebRecHit.et,"et[n]/F");
436			ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
437			ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
438	yilmaz	1.4	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
439	yilmaz	1.2
440	yilmaz	1.7	towerTree = fs->make<TTree>("tower","");
441			towerTree->Branch("n",&myTowers.n,"n/I");
442			towerTree->Branch("e",myTowers.e,"e[n]/F");
443			towerTree->Branch("et",myTowers.et,"et[n]/F");
444			towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
445			towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
446			towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
447
448
449	yilmaz	1.6	if(doBasicClusters_){
450			bcTree = fs->make<TTree>("bc","");
451			bcTree->Branch("n",&myBC.n,"n/I");
452			bcTree->Branch("e",myBC.e,"e[n]/F");
453			bcTree->Branch("et",myBC.et,"et[n]/F");
454			bcTree->Branch("eta",myBC.eta,"eta[n]/F");
455			bcTree->Branch("phi",myBC.phi,"phi[n]/F");
456			bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
457			bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
458			bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
459			// bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
460			}
461
462	yilmaz	1.7	if(doFastJet_){
463			bkgTree = fs->make<TTree>("bkg","");
464			bkgTree->Branch("n",&bkg.n,"n/I");
465			bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
466			bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
467			}
468	yilmaz	1.6
469	yilmaz	1.1	}
470
471			// ------------ method called once each job just after ending the event loop ------------
472			void
473			RecHitTreeProducer::endJob() {
474			}
475
476			double RecHitTreeProducer::getEt(const DetId &id, double energy){
477			const GlobalPoint& pos=geo->getPosition(id);
478			double et = energy*sin(pos.theta());
479			return et;
480			}
481	nart	1.3
482	yilmaz	1.1	double RecHitTreeProducer::getEta(const DetId &id){
483			const GlobalPoint& pos=geo->getPosition(id);
484			double et = pos.eta();
485			return et;
486			}
487
488			double RecHitTreeProducer::getPhi(const DetId &id){
489			const GlobalPoint& pos=geo->getPosition(id);
490			double et = pos.phi();
491			return et;
492			}
493
494
495
496			//define this as a plug-in
497			DEFINE_FWK_MODULE(RecHitTreeProducer);