// -*- C++ -*-
//
// Package:    HFJetsAnalyzer
// Class:      HFJetsAnalyzer
// 
// class HFJetsAnalyzer HFJetsAnalyzer.cc analyzer/HFJetsAnalyzer/src/HFJetsAnalyzer.cc
// Original Author:  Andrey Pozdnyakov
//         Created:  Mon Sep 24 13:31:07 CEST 2007

// system include files

#include <memory>

#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/ParameterSet/interface/ParameterSet.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "RecoJets/JetAlgorithms/interface/JetAlgoHelper.h"
#include "DataFormats/JetReco/interface/CaloJetCollection.h"
#include "DataFormats/JetReco/interface/CaloJet.h"
#include "DataFormats/JetReco/interface/GenJet.h"
#include <TROOT.h>
#include <TFile.h>
#include <TH1.h>
#include <TProfile.h>
#include <cmath>
using namespace edm;
using namespace reco;
using namespace std;

class HFJetsAnalyzer : public edm::EDAnalyzer {
   public:
      explicit HFJetsAnalyzer(const edm::ParameterSet&);
      ~HFJetsAnalyzer();

   private:
      virtual void beginJob(const edm::EventSetup&) ;
      virtual void analyze(const edm::Event&, const edm::EventSetup&);
      virtual void endJob() ;

  string CaloJetAlgorithm, GenJetAlgorithm;
  string outputFileName_;
  double etamaxCut_, etaminCut_, ptmaxCut_, ptminCut_, coneSize1_, coneSize2_;
  TH1F h_numGen, h_numCal;
  TH1F h_ptCal, h_etaCal, h_phiCal;
  TH1F h_ptGen, h_etaGen, h_phiGen;
  TH1F h_dRminGenCalo, h_numCaloInCone;
  TH1F h_dRminCaloGen, h_numGenInCone;
  TProfile  p_CaloEnRespVsEta,  p_CaloPtRespVsEta;
  TProfile  p_CaloEnRespVsEta2,  p_CaloPtRespVsEta2;
  TFile* m_file;
};

HFJetsAnalyzer::HFJetsAnalyzer(const edm::ParameterSet& cfg)
{
  CaloJetAlgorithm = cfg.getParameter<string>("CaloJetAlgorithm");
  GenJetAlgorithm  = cfg.getParameter<string>("GenJetAlgorithm");
  etamaxCut_ = cfg.getParameter<double>("etamaxCut");
  etaminCut_ = cfg.getParameter<double>("etaminCut");
  ptmaxCut_ = cfg.getParameter<double>("ptmaxCut");
  ptminCut_ = cfg.getParameter<double>("ptminCut");
  coneSize1_ = cfg.getParameter<double>("coneSize1");
  coneSize2_ = cfg.getParameter<double>("coneSize2");
  outputFileName_ = cfg.getParameter<string>("outputFileName");
}

HFJetsAnalyzer::~HFJetsAnalyzer()
{}

void HFJetsAnalyzer::beginJob(const edm::EventSetup&)
{
// Open the histogram file and book some associated histograms
  m_file   = new TFile(outputFileName_.c_str(),"RECREATE");
  h_numGen = TH1F("numGen", "Number of HF GenJets in event", 30, 0, 30);
  h_numCal = TH1F("numCal", "Number of HF CaloJets in event", 30, 0, 30);
  h_ptCal  = TH1F("ptCal",  "p_{T} of leading CaloJets", 10, 0, ptmaxCut_);
  h_etaCal = TH1F("etaCal", "#eta of CaloJets", 50, -6, 6 );
  h_phiCal = TH1F("phiCal", "#phi of CaloJets", 50, -M_PI, M_PI);
  h_ptGen  = TH1F("ptGen",  "p_{T} of GenJets", 10, 0, ptmaxCut_);
  h_etaGen = TH1F("etaGen", "#eta of GenJets", 50, -6, 6 );
  h_phiGen = TH1F("phiGen", "#phi of GenJets", 50, -M_PI, M_PI);
  h_dRminGenCalo = TH1F("dRminGenCalo","Min #DeltaR from Gen to Calo jet", 50, 0, 2);
  h_numCaloInCone = TH1F("numCaloInCone","Number of Calo jets in coneSize1 near Gen jet", 5, 0, 5);
  p_CaloEnRespVsEta = TProfile("CaloEnRespVsEta","Jets Energy Response vs Eta", 50, etaminCut_, etamaxCut_);
  p_CaloPtRespVsEta = TProfile("CaloPtRespVsEta","Jets Pt Response vs Eta", 50, etaminCut_, etamaxCut_);

  h_dRminCaloGen = TH1F("dRminCaloGen","Min #DeltaR from Cal to Gen jet", 50, 0, 2);
  h_numGenInCone = TH1F("numGenInCone","Number of Gen jets in coneSize2 near Calo jet", 5, 0, 5);
  p_CaloEnRespVsEta2 = TProfile("CaloEnRespVsEta2","Jets Energy Response vs Eta", 50, etaminCut_, etamaxCut_);
  p_CaloPtRespVsEta2 = TProfile("CaloPtRespVsEta2","Jets Pt Response vs Eta", 50, etaminCut_, etamaxCut_);
}

void HFJetsAnalyzer::analyze(const edm::Event& evt, const edm::EventSetup& es)
{
  Handle<CaloJetCollection> caloJets;
  evt.getByLabel(CaloJetAlgorithm, caloJets);

//First, loop over all Calo jets (appropriate to cuts) and fill histogramms
  int jetInd = 0;
  for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++cal) 
    {
      if (abs(cal->eta()) > etaminCut_ && abs(cal->eta()) < etamaxCut_ && cal->pt() > ptminCut_ && cal->pt() < ptmaxCut_)
	{
	  //std::cout << "CALO JET #" << jetInd << std::endl << cal->print() << std::endl;
	  h_ptCal.Fill( cal->pt() );
	  h_etaCal.Fill( cal->eta() );
	  h_phiCal.Fill( cal->phi() );
	  jetInd++;
	}
    }   
  h_numCal.Fill(jetInd);

  Handle<GenJetCollection> genJets;
  evt.getByLabel(GenJetAlgorithm, genJets);

//Second, loop over Gen jets (with same cuts) and fill same histograms as for Calo jets
  jetInd = 0;
  for( GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end(); ++gen)
    {
      if (abs(gen->eta()) > etaminCut_ && abs(gen->eta()) < etamaxCut_ && gen->pt() > ptminCut_ && gen->pt() < ptmaxCut_) 
	{ 
	  //std::cout << "GEN JET #" << jetInd << std::endl << gen->print() << std::endl;
	  h_ptGen.Fill( gen->pt() );
	  h_etaGen.Fill( gen->eta() );
	  h_phiGen.Fill( gen->phi() );
	  jetInd++;
	  
	}
    }
  h_numGen.Fill(jetInd);

//Now loop over Gen jets and find nearest Calo jets in "coneSize1"
  jetInd = 0;
  for( GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end(); ++gen)
    {
      if (abs(gen->eta()) > etaminCut_ && abs(gen->eta()) < etamaxCut_ && gen->pt() > ptminCut_ && gen->pt() < ptmaxCut_) 
	{ 
	  //std::cout << "GEN JET #" << jetInd << std::endl << gen->print() << std::endl;
	  h_ptGen.Fill( gen->pt() );
	  h_etaGen.Fill( gen->eta() );
	  h_phiGen.Fill( gen->phi() );
	  jetInd++;
	  
	  CaloJet MatchedJet;	
	  
	  double dRminSave=1000.0;
	  int JetsInCone=0;
	  for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end() && cal->pt()>ptminCut_; ++ cal ) 
	    { 
	      double delR = deltaR( cal->eta(), cal->phi(), gen->eta(), gen->phi() ); 
	      if(delR<dRminSave)
		{
		  dRminSave=delR;
		  MatchedJet = *cal;
		}
	      h_dRminGenCalo.Fill(dRminSave);
	      if(delR < coneSize1_)
		{
		  JetsInCone++;
		}
	    }
	  p_CaloEnRespVsEta.Fill(fabs(MatchedJet.eta()),MatchedJet.energy()/gen->energy());		
	  p_CaloPtRespVsEta.Fill(fabs(MatchedJet.eta()),MatchedJet.pt()/gen->pt());		
	  h_numCaloInCone.Fill(JetsInCone);
	}
    }
//Now do the opposite thing: loop over Cal jets and find nearest Gen jets in "coneSize2"
  jetInd = 0;
  for( CaloJetCollection::const_iterator cal = caloJets->begin(); cal != caloJets->end(); ++cal)
    {
      if (abs(cal->eta()) > etaminCut_ && abs(cal->eta()) < etamaxCut_ && cal->pt() > ptminCut_ && cal->pt() < ptmaxCut_) 
	{ 
	  //std::cout << "GEN JET #" << jetInd << std::endl << gen->print() << std::endl;
	  h_ptGen.Fill( cal->pt() );
	  h_etaGen.Fill( cal->eta() );
	  h_phiGen.Fill( cal->phi() );
	  jetInd++;
	  
	  GenJet MatchedJet;	
	  
	  double dRminSave=1000.0;
	  int JetsInCone=0;
	  for( GenJetCollection::const_iterator gen = genJets->begin(); gen != genJets->end() && gen->pt()>ptminCut_; ++ gen) 
	    { 
	      double delR = deltaR( cal->eta(), cal->phi(), gen->eta(), gen->phi() ); 
	      if(delR<dRminSave)
		{
		  dRminSave=delR;
		  MatchedJet = *gen;
		}
	      h_dRminCaloGen.Fill(dRminSave);
	      if(delR < coneSize2_) 
		{
		  JetsInCone++;
		}
	    }
	  p_CaloEnRespVsEta2.Fill(fabs(MatchedJet.eta()), cal->energy()/MatchedJet.energy());		
	  p_CaloPtRespVsEta2.Fill(fabs(MatchedJet.eta()), cal->pt()/MatchedJet.pt());		
	  h_numGenInCone.Fill(JetsInCone);
	}
    }
}
void HFJetsAnalyzer::endJob() 
{
  m_file ->Write();
}

//define this as a plug-in
DEFINE_FWK_MODULE(HFJetsAnalyzer);