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.9 2010/11/09 22:54:45 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 depth[MAXHITS];
  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;

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

  double cone;
  double hfTowerThreshold_;
  double hfLongThreshold_;
  double hfShortThreshold_;
  double hbheThreshold_;
  double ebThreshold_;
  double eeThreshold_;

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

  bool doEbyEonly_;
};

//
// 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("iterativeConePu5CaloJets"));
  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"));
  doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
  hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
  hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",3.);
  hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",3.);
}


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

   int nHFlong = 0;
   int nHFshort = 0;
   int nHFtower = 0;
   
   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;
     hfRecHit.depth[hfRecHit.n] = hit.id().depth();
     if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshort++;
     if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlong++;

     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++;
   }
   if(!doEbyEonly_){
   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_ && !doEbyEonly_){
   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(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtower++;

      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_ && !doEbyEonly_){
      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(); 
      }
   }

   if(!doEbyEonly_){
     towerTree->Fill();
     
     eeTree->Fill();
     ebTree->Fill();
     
     hbheTree->Fill();
     hfTree->Fill();
      
     if (doFastJet_) {
       bkgTree->Fill();
     }
   }

   nt->Fill(nHFtower,nHFlong,nHFshort);
   
}


// ------------ 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("depth",hfRecHit.depth,"depth[n]/I");
  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");
  }
  
  nt = fs->make<TNtuple>("ntEvent","","nHF:nHFlong:nHFshort");

}

// ------------ 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.10
Committed:	Tue Nov 9 22:55:50 2010 UTC (14 years, 5 months ago) by yilmaz
Content type:	text/plain
Branch:	MAIN
Changes since 1.9:	+2 -2 lines
Log Message:	event by event ntuple
#	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.10	// $Id: RecHitTreeProducer.cc,v 1.9 2010/11/09 22:54:45 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	yilmaz	1.9	int depth[MAXHITS];
72	yilmaz	1.1	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.9	TNtuple* nt;
144	yilmaz	1.1	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.9	double hfTowerThreshold_;
154			double hfLongThreshold_;
155			double hfShortThreshold_;
156			double hbheThreshold_;
157			double ebThreshold_;
158			double eeThreshold_;
159	yilmaz	1.1
160			edm::Service<TFileService> fs;
161			const CentralityBins * cbins_;
162			const CaloGeometry *geo;
163
164			edm::InputTag HcalRecHitHFSrc_;
165			edm::InputTag HcalRecHitHBHESrc_;
166	yilmaz	1.2	edm::InputTag EBSrc_;
167			edm::InputTag EESrc_;
168	yilmaz	1.6	edm::InputTag BCSrc_;
169			edm::InputTag TowerSrc_;
170
171	nart	1.3	edm::InputTag JetSrc_;
172	yilmaz	1.6
173	yilmaz	1.7	edm::InputTag FastJetTag_;
174
175	yilmaz	1.6	bool useJets_;
176			bool doBasicClusters_;
177	yilmaz	1.7	bool doTowers_;
178			bool doEcal_;
179			bool doHcal_;
180	yilmaz	1.1
181	yilmaz	1.7	bool doFastJet_;
182	yilmaz	1.9
183			bool doEbyEonly_;
184	yilmaz	1.1	};
185
186			//
187			// constants, enums and typedefs
188			//
189
190			//
191			// static data member definitions
192			//
193
194			//
195			// constructors and destructor
196			//
197			RecHitTreeProducer::RecHitTreeProducer(const edm::ParameterSet& iConfig) :
198			cbins_(0),
199			geo(0),
200			cone(0.5)
201			{
202			//now do what ever initialization is needed
203	yilmaz	1.2	EBSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EBRecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEB"));
204			EESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("EERecHitSrc",edm::InputTag("ecalRecHit","EcalRecHitsEE"));
205	yilmaz	1.1	HcalRecHitHFSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHFRecHitSrc",edm::InputTag("hfreco"));
206			HcalRecHitHBHESrc_ = iConfig.getUntrackedParameter<edm::InputTag>("hcalHBHERecHitSrc",edm::InputTag("hbhereco"));
207	yilmaz	1.6	BCSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("BasicClusterSrc1",edm::InputTag("ecalRecHit","EcalRecHitsEB","RECO"));
208			TowerSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("towersSrc",edm::InputTag("towerMaker"));
209	yilmaz	1.10	JetSrc_ = iConfig.getUntrackedParameter<edm::InputTag>("JetSrc",edm::InputTag("iterativeConePu5CaloJets"));
210	yilmaz	1.7	useJets_ = iConfig.getUntrackedParameter<bool>("useJets",true);
211	yilmaz	1.6	doBasicClusters_ = iConfig.getUntrackedParameter<bool>("doBasicClusters",false);
212	yilmaz	1.7	doTowers_ = iConfig.getUntrackedParameter<bool>("doTowers",true);
213			doEcal_ = iConfig.getUntrackedParameter<bool>("doEcal",true);
214			doHcal_ = iConfig.getUntrackedParameter<bool>("doHcal",true);
215			doFastJet_ = iConfig.getUntrackedParameter<bool>("doFastJet",true);
216			FastJetTag_ = iConfig.getUntrackedParameter<edm::InputTag>("FastJetTag",edm::InputTag("kt4CaloJets"));
217	yilmaz	1.9	doEbyEonly_ = iConfig.getUntrackedParameter<bool>("doEbyEonly",false);
218			hfTowerThreshold_ = iConfig.getUntrackedParameter<double>("HFtowerMin",3.);
219			hfLongThreshold_ = iConfig.getUntrackedParameter<double>("HFlongMin",3.);
220			hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",3.);
221	yilmaz	1.1	}
222
223
224			RecHitTreeProducer::~RecHitTreeProducer()
225			{
226
227			// do anything here that needs to be done at desctruction time
228			// (e.g. close files, deallocate resources etc.)
229
230			}
231
232
233			//
234			// member functions
235			//
236
237			// ------------ method called to for each event ------------
238			void
239			RecHitTreeProducer::analyze(const edm::Event& ev, const edm::EventSetup& iSetup)
240			{
241
242			hfRecHit.n = 0;
243			hbheRecHit.n = 0;
244	yilmaz	1.2	ebRecHit.n = 0;
245			eeRecHit.n = 0;
246	yilmaz	1.7	myBC.n = 0;
247			myTowers.n = 0;
248			bkg.n = 0;
249	yilmaz	1.1
250	yilmaz	1.7	if(doEcal_){
251	yilmaz	1.2	ev.getByLabel(EBSrc_,ebHits);
252			ev.getByLabel(EESrc_,eeHits);
253	yilmaz	1.7	}
254			if(doHcal_){
255	yilmaz	1.1	ev.getByLabel(HcalRecHitHFSrc_,hfHits);
256			ev.getByLabel(HcalRecHitHBHESrc_,hbheHits);
257	yilmaz	1.7	}
258	yilmaz	1.6	if(useJets_) {
259			ev.getByLabel(JetSrc_,jets);
260			}
261	nart	1.3
262	yilmaz	1.6	if(doBasicClusters_){
263			ev.getByLabel(BCSrc_,bClusters);
264	nart	1.3	}
265	yilmaz	1.6
266	yilmaz	1.7	if(doTowers_){
267			ev.getByLabel(TowerSrc_,towers);
268			}
269	yilmaz	1.6
270	yilmaz	1.7	if(doFastJet_){
271			ev.getByLabel(edm::InputTag(FastJetTag_.label(),"rhos",FastJetTag_.process()),rhos);
272			ev.getByLabel(edm::InputTag(FastJetTag_.label(),"sigmas",FastJetTag_.process()),sigmas);
273			bkg.n = rhos->size();
274			for(unsigned int i = 0; i < rhos->size(); ++i){
275			bkg.rho[i] = (*rhos)[i];
276			bkg.sigma[i] = (*sigmas)[i];
277			}
278			}
279	nart	1.3
280	yilmaz	1.6	if(0 && !cbins_) cbins_ = getCentralityBinsFromDB(iSetup);
281	yilmaz	1.1
282			if(!geo){
283			edm::ESHandle<CaloGeometry> pGeo;
284			iSetup.get<CaloGeometryRecord>().get(pGeo);
285			geo = pGeo.product();
286			}
287	yilmaz	1.9
288			int nHFlong = 0;
289			int nHFshort = 0;
290			int nHFtower = 0;
291	nart	1.3
292	yilmaz	1.7	if(doHcal_){
293	yilmaz	1.1	for(unsigned int i = 0; i < hfHits->size(); ++i){
294			const HFRecHit & hit= (*hfHits)[i];
295			hfRecHit.e[hfRecHit.n] = hit.energy();
296			hfRecHit.et[hfRecHit.n] = getEt(hit.id(),hit.energy());
297			hfRecHit.eta[hfRecHit.n] = getEta(hit.id());
298			hfRecHit.phi[hfRecHit.n] = getPhi(hit.id());
299	nart	1.3	hfRecHit.isjet[hfRecHit.n] = false;
300	yilmaz	1.9	hfRecHit.depth[hfRecHit.n] = hit.id().depth();
301			if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshort++;
302			if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlong++;
303
304	yilmaz	1.7	if(useJets_){
305	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
306			const reco::Jet& jet = (*jets)[j];
307			double dr = reco::deltaR(hfRecHit.eta[hfRecHit.n],hfRecHit.phi[hfRecHit.n],jet.eta(),jet.phi());
308			if(dr < cone){ hfRecHit.isjet[hfRecHit.n] = true; }
309			}
310			}
311	yilmaz	1.1	hfRecHit.n++;
312			}
313	yilmaz	1.9	if(!doEbyEonly_){
314	yilmaz	1.1	for(unsigned int i = 0; i < hbheHits->size(); ++i){
315			const HBHERecHit & hit= (*hbheHits)[i];
316	nart	1.3	hbheRecHit.e[hbheRecHit.n] = hit.energy();
317			hbheRecHit.et[hbheRecHit.n] = getEt(hit.id(),hit.energy());
318			hbheRecHit.eta[hbheRecHit.n] = getEta(hit.id());
319			hbheRecHit.phi[hbheRecHit.n] = getPhi(hit.id());
320			hbheRecHit.isjet[hbheRecHit.n] = false;
321	yilmaz	1.7	if(useJets_){
322	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
323			const reco::Jet& jet = (*jets)[j];
324			double dr = reco::deltaR(hbheRecHit.eta[hbheRecHit.n],hbheRecHit.phi[hbheRecHit.n],jet.eta(),jet.phi());
325			if(dr < cone){ hbheRecHit.isjet[hbheRecHit.n] = true; }
326			}
327			}
328	yilmaz	1.2	hbheRecHit.n++;
329			}
330	yilmaz	1.7	}
331	yilmaz	1.9	}
332			if(doEcal_ && !doEbyEonly_){
333	yilmaz	1.2	for(unsigned int i = 0; i < ebHits->size(); ++i){
334			const EcalRecHit & hit= (*ebHits)[i];
335			ebRecHit.e[ebRecHit.n] = hit.energy();
336			ebRecHit.et[ebRecHit.n] = getEt(hit.id(),hit.energy());
337			ebRecHit.eta[ebRecHit.n] = getEta(hit.id());
338			ebRecHit.phi[ebRecHit.n] = getPhi(hit.id());
339	nart	1.3	ebRecHit.isjet[ebRecHit.n] = false;
340	yilmaz	1.7	if(useJets_){
341	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
342			const reco::Jet& jet = (*jets)[j];
343			double dr = reco::deltaR(ebRecHit.eta[ebRecHit.n],ebRecHit.phi[ebRecHit.n],jet.eta(),jet.phi());
344			if(dr < cone){ ebRecHit.isjet[ebRecHit.n] = true; }
345			}
346			}
347	yilmaz	1.2	ebRecHit.n++;
348			}
349	nart	1.3
350	yilmaz	1.2	for(unsigned int i = 0; i < eeHits->size(); ++i){
351			const EcalRecHit & hit= (*eeHits)[i];
352			eeRecHit.e[eeRecHit.n] = hit.energy();
353			eeRecHit.et[eeRecHit.n] = getEt(hit.id(),hit.energy());
354			eeRecHit.eta[eeRecHit.n] = getEta(hit.id());
355			eeRecHit.phi[eeRecHit.n] = getPhi(hit.id());
356	nart	1.3	eeRecHit.isjet[eeRecHit.n] = false;
357	yilmaz	1.7	if(useJets_){
358	nart	1.3	for(unsigned int j = 0 ; j < jets->size(); ++j){
359			const reco::Jet& jet = (*jets)[j];
360			double dr = reco::deltaR(eeRecHit.eta[eeRecHit.n],eeRecHit.phi[eeRecHit.n],jet.eta(),jet.phi());
361			if(dr < cone){ eeRecHit.isjet[eeRecHit.n] = true; }
362			}
363			}
364	yilmaz	1.2	eeRecHit.n++;
365	yilmaz	1.1	}
366	yilmaz	1.7	}
367
368			if(doTowers_){
369
370			for(unsigned int i = 0; i < towers->size(); ++i){
371			const CaloTower & hit= (*towers)[i];
372			myTowers.e[myTowers.n] = hit.energy();
373			myTowers.et[myTowers.n] = getEt(hit.id(),hit.energy());
374			myTowers.eta[myTowers.n] = getEta(hit.id());
375			myTowers.phi[myTowers.n] = getPhi(hit.id());
376			myTowers.isjet[myTowers.n] = false;
377	yilmaz	1.9
378			if(hit.ieta() > 29 && hit.energy() > hfTowerThreshold_) nHFtower++;
379
380	yilmaz	1.7	if(useJets_){
381			for(unsigned int j = 0 ; j < jets->size(); ++j){
382			const reco::Jet& jet = (*jets)[j];
383			double dr = reco::deltaR(myTowers.eta[myTowers.n],myTowers.phi[myTowers.n],jet.eta(),jet.phi());
384			if(dr < cone){ myTowers.isjet[myTowers.n] = true; }
385			}
386			}
387			myTowers.n++;
388			}
389	yilmaz	1.6
390	yilmaz	1.7	}
391	yilmaz	1.6
392	yilmaz	1.9	if(doBasicClusters_ && !doEbyEonly_){
393	yilmaz	1.6	for(unsigned int j = 0 ; j < jets->size(); ++j){
394			const reco::Jet& jet = (*jets)[j];
395			myBC.n = 0;
396			myBC.jtpt = jet.pt();
397			myBC.jteta = jet.eta();
398			myBC.jtphi = jet.phi();
399
400			for(unsigned int i = 0; i < bClusters->size(); ++i){
401			const reco::BasicCluster & bc= (*bClusters)[i];
402			double dr = reco::deltaR(bc.eta(),bc.phi(),jet.eta(),jet.phi());
403			if(dr < cone){
404			myBC.e[myBC.n] = bc.energy();
405			myBC.et[myBC.n] = bc.energy()*sin(bc.position().theta());
406			myBC.eta[myBC.n] = bc.eta();
407			myBC.phi[myBC.n] = bc.phi();
408			myBC.n++;
409			}
410			}
411			bcTree->Fill();
412			}
413			}
414	yilmaz	1.7
415	yilmaz	1.9	if(!doEbyEonly_){
416			towerTree->Fill();
417
418			eeTree->Fill();
419			ebTree->Fill();
420
421			hbheTree->Fill();
422			hfTree->Fill();
423
424			if (doFastJet_) {
425			bkgTree->Fill();
426			}
427			}
428
429			nt->Fill(nHFtower,nHFlong,nHFshort);
430	nart	1.3
431	yilmaz	1.1	}
432
433
434			// ------------ method called once each job just before starting event loop ------------
435			void
436			RecHitTreeProducer::beginJob()
437			{
438	nart	1.3
439	yilmaz	1.1	hbheTree = fs->make<TTree>("hbhe","");
440			hbheTree->Branch("n",&hbheRecHit.n,"n/I");
441			hbheTree->Branch("e",hbheRecHit.e,"e[n]/F");
442			hbheTree->Branch("et",hbheRecHit.et,"et[n]/F");
443			hbheTree->Branch("eta",hbheRecHit.eta,"eta[n]/F");
444			hbheTree->Branch("phi",hbheRecHit.phi,"phi[n]/F");
445	yilmaz	1.5	hbheTree->Branch("isjet",hbheRecHit.isjet,"isjet[n]/O");
446	nart	1.3
447	yilmaz	1.1	hfTree = fs->make<TTree>("hf","");
448			hfTree->Branch("n",&hfRecHit.n,"n/I");
449			hfTree->Branch("e",hfRecHit.e,"e[n]/F");
450			hfTree->Branch("et",hfRecHit.et,"et[n]/F");
451			hfTree->Branch("eta",hfRecHit.eta,"eta[n]/F");
452			hfTree->Branch("phi",hfRecHit.phi,"phi[n]/F");
453	yilmaz	1.9	hfTree->Branch("depth",hfRecHit.depth,"depth[n]/I");
454	yilmaz	1.5	hfTree->Branch("isjet",hfRecHit.isjet,"isjet[n]/O");
455	yilmaz	1.1
456	yilmaz	1.2	eeTree = fs->make<TTree>("ee","");
457			eeTree->Branch("n",&eeRecHit.n,"n/I");
458			eeTree->Branch("e",eeRecHit.e,"e[n]/F");
459			eeTree->Branch("et",eeRecHit.et,"et[n]/F");
460			eeTree->Branch("eta",eeRecHit.eta,"eta[n]/F");
461			eeTree->Branch("phi",eeRecHit.phi,"phi[n]/F");
462	yilmaz	1.5	eeTree->Branch("isjet",eeRecHit.isjet,"isjet[n]/O");
463	nart	1.3
464	yilmaz	1.2	ebTree = fs->make<TTree>("eb","");
465			ebTree->Branch("n",&ebRecHit.n,"n/I");
466			ebTree->Branch("e",ebRecHit.e,"e[n]/F");
467			ebTree->Branch("et",ebRecHit.et,"et[n]/F");
468			ebTree->Branch("eta",ebRecHit.eta,"eta[n]/F");
469			ebTree->Branch("phi",ebRecHit.phi,"phi[n]/F");
470	yilmaz	1.4	ebTree->Branch("isjet",ebRecHit.isjet,"isjet[n]/O");
471	yilmaz	1.2
472	yilmaz	1.7	towerTree = fs->make<TTree>("tower","");
473			towerTree->Branch("n",&myTowers.n,"n/I");
474			towerTree->Branch("e",myTowers.e,"e[n]/F");
475			towerTree->Branch("et",myTowers.et,"et[n]/F");
476			towerTree->Branch("eta",myTowers.eta,"eta[n]/F");
477			towerTree->Branch("phi",myTowers.phi,"phi[n]/F");
478			towerTree->Branch("isjet",myTowers.isjet,"isjet[n]/O");
479
480
481	yilmaz	1.6	if(doBasicClusters_){
482			bcTree = fs->make<TTree>("bc","");
483			bcTree->Branch("n",&myBC.n,"n/I");
484			bcTree->Branch("e",myBC.e,"e[n]/F");
485			bcTree->Branch("et",myBC.et,"et[n]/F");
486			bcTree->Branch("eta",myBC.eta,"eta[n]/F");
487			bcTree->Branch("phi",myBC.phi,"phi[n]/F");
488			bcTree->Branch("jtpt",&myBC.jtpt,"jtpt/F");
489			bcTree->Branch("jteta",&myBC.jteta,"jteta/F");
490			bcTree->Branch("jtphi",&myBC.jtphi,"jtphi/F");
491			// bcTree->Branch("isjet",bcRecHit.isjet,"isjet[n]/O");
492			}
493
494	yilmaz	1.7	if(doFastJet_){
495			bkgTree = fs->make<TTree>("bkg","");
496			bkgTree->Branch("n",&bkg.n,"n/I");
497			bkgTree->Branch("rho",bkg.rho,"rho[n]/F");
498			bkgTree->Branch("sigma",bkg.sigma,"sigma[n]/F");
499			}
500	yilmaz	1.9
501			nt = fs->make<TNtuple>("ntEvent","","nHF:nHFlong:nHFshort");
502	yilmaz	1.6
503	yilmaz	1.1	}
504
505			// ------------ method called once each job just after ending the event loop ------------
506			void
507			RecHitTreeProducer::endJob() {
508			}
509
510			double RecHitTreeProducer::getEt(const DetId &id, double energy){
511			const GlobalPoint& pos=geo->getPosition(id);
512			double et = energy*sin(pos.theta());
513			return et;
514			}
515	nart	1.3
516	yilmaz	1.1	double RecHitTreeProducer::getEta(const DetId &id){
517			const GlobalPoint& pos=geo->getPosition(id);
518			double et = pos.eta();
519			return et;
520			}
521
522			double RecHitTreeProducer::getPhi(const DetId &id){
523			const GlobalPoint& pos=geo->getPosition(id);
524			double et = pos.phi();
525			return et;
526			}
527
528
529
530			//define this as a plug-in
531			DEFINE_FWK_MODULE(RecHitTreeProducer);