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.11 2010/11/09 22:57: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 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",0.5);
  hfShortThreshold_ = iConfig.getUntrackedParameter<double>("HFshortMin",0.85);
}


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 nHFlongPlus = 0;
   int nHFshortPlus = 0;
   int nHFtowerPlus = 0;
   int nHFlongMinus = 0;
   int nHFshortMinus = 0;
   int nHFtowerMinus = 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.id().ieta() > 0){
     if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortPlus++;
     if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongPlus++;
     }else{
       if(hit.energy() > hfShortThreshold_ && hit.id().depth() != 1) nHFshortMinus++;
       if(hit.energy() > hfLongThreshold_ && hit.id().depth() == 1) nHFlongMinus++;
     }

     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_) nHFtowerPlus++;
      if(hit.ieta() < -29 && hit.energy() > hfTowerThreshold_) nHFtowerMinus++;

      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(nHFtowerPlus,nHFtowerMinus,nHFlongPlus,nHFlongMinus,nHFshortPlus,nHFshortMinus);
   
}


// ------------ 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","","nHFplus:nHFminus:nHFlongPlus:nHFlongMinus:nHFshortPlus:nHFshortMinus");

}

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