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.12 2010/11/11 18:37:16 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_;
  
  double hbhePtMin_;
  double hfPtMin_;
  double ebPtMin_;
  double eePtMin_;
  double towerPtMin_;
  
   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);
  hbhePtMin_ = iConfig.getUntrackedParameter<double>("HBHETreePtMin",0);
  hfPtMin_ = iConfig.getUntrackedParameter<double>("HFTreePtMin",0);
  ebPtMin_ = iConfig.getUntrackedParameter<double>("EBTreePtMin",0);
  eePtMin_ = iConfig.getUntrackedParameter<double>("EETreePtMin",0.);
  towerPtMin_ = iConfig.getUntrackedParameter<double>("TowerTreePtMin",0.);
}


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