#include "Rivet/Analysis.hh"
#include "Rivet/Run.hh"
#include "HepMC/GenEvent.h"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/PromptFinalState.hh"
#include "Rivet/Projections/UnstableFinalState.hh"
#include "Rivet/Projections/VetoedFinalState.hh"
#include "Rivet/Projections/ChargedLeptons.hh"
#include "Rivet/Projections/DressedLeptons.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/AnalysisLoader.hh"
#include "fastjet/tools/Filter.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/AreaDefinition.hh"
#include "fastjet/tools/JetMedianBackgroundEstimator.hh"
#include "fastjet/tools/Subtractor.hh"
//fastjet::GhostedAreaSpec areaspec(2.5,1,0.01);
//fastjet::AreaDefinition area_def(fastjet::active_area_explicit_ghosts,areaspec);
namespace Rivet {

  class tt_signal : public Analysis {
  public:

    /// Minimal constructor
    tt_signal() : Analysis("tt_signal"), areaspec(new fastjet::GhostedAreaSpec(ymax, nrepeat, ghost_area)), area_def(new fastjet::AreaDefinition(fastjet::active_area_explicit_ghosts, *areaspec))
    {
    }
    ~tt_signal()
    {
     delete area_def; 
     delete areaspec; 
    }
/***************************************************************************************/
// Compute neutrino attributes for reconstructing the leptonic top
    double computeneutrinoz(const FourMomentum& lepton, FourMomentum& met) const {
      //computing z component of neutrino momentum given lepton and met
      double pzneutrino;
      double m_W = 80.399; // in GeV, given in the paper
      double k = (( sqr( m_W ) - sqr( lepton.mass() ) ) / 2 ) + (lepton.px() * met.px() + lepton.py() * met.py());
      double a = sqr ( lepton.E() )- sqr ( lepton.pz() );
      double b = -2*k*lepton.pz();
      double c = sqr( lepton.E() ) * sqr( met.pT() ) - sqr( k );
      double discriminant = sqr(b) - 4 * a * c;
      double quad[2] = { (- b - sqrt(discriminant)) / (2 * a), (- b + sqrt(discriminant)) / (2 * a) }; //two possible quadratic solns
      if (discriminant < 0)  pzneutrino = - b / (2 * a); //if the discriminant is negative
      else { //if the discriminant is greater than or equal to zero, take the soln with smallest absolute value
        double absquad[2];
        for (int n=0; n<2; ++n)  absquad[n] = fabs(quad[n]);
        if (absquad[0] < absquad[1])  pzneutrino = quad[0];
        else                          pzneutrino = quad[1];
      }
      if ( !std::isfinite(pzneutrino) )  std::cout << "Found non-finite value" << std::endl;
      return pzneutrino;
    }
/***************************************************************************************/
// Reconstruct top quark masses and pT and then reconstruct Z' mass
    void PartonLevel(PseudoJets st_jets, const Event& event, int ind) {
      const double weight = event.weight();
     
     // RT Trimming procedure (take out small jets out of the large jet if pT < fcut*pT(large)
     if (ind==0) {
      for (int i=0; i < st_jets.size(); i++) {
        PseudoJet copy_jets;
	vector<PseudoJet> c_constits;
        PseudoJets constits = fastjet::sorted_by_pt(st_jets[i].constituents());
        foreach (const PseudoJet& cjet, constits) {
	  if (cjet.pt() > fcut*st_jets[i].pt()) c_constits.push_back(cjet); //copy += cjet;
        }
	copy_jets = fastjet::join(c_constits);
	st_jets[i] = copy_jets;
      }
     }
     // cutoffs
      PseudoJets trimmed_fat_jets;
      for (size_t i = 0; i < st_jets.size(); ++i) {
        const Jet tj = st_jets[i];
        if (tj.mass()          <= mass_cutoff*GeV) continue;
        if (tj.pt()            <= pT_cutoff*GeV) continue;
        if (tj.abseta()        >=     2.0) continue;
        trimmed_fat_jets.push_back(tj);
      }
      if (trimmed_fat_jets.empty()) {
	numVetoedEvents[ind] += 1;
	numEmpty[ind] += 1;
	vetoEvent;
      }
      st_jets = trimmed_fat_jets;

      // Bare leptons 
      const Particles& bare_lep = apply<PromptFinalState>(event, "bare_lep").particles();
      const Particles& bare_tau = apply<PromptFinalState>(event, "bare_tau").particles();
      if (bare_lep.size() + bare_tau.size() != 1) {
	numVetoedEvents[ind] += 1;
	numEarlyReject[ind] += 1;
	vetoEvent;
      }
      // Electrons and muons
      const vector<DressedLepton>& electrons = apply<DressedLeptons>(event, "electrons").dressedLeptons();
      const vector<DressedLepton>& muons = apply<DressedLeptons>(event, "muons").dressedLeptons();
      if (electrons.size() + muons.size() != 1) {
        numVetoedEvents[ind] += 1;
	numEarlyReject[ind] += 1;
	vetoEvent;
	}
      const DressedLepton& lepton = muons.empty() ? electrons[0] : muons[0];
      // Get the neutrinos from the event record (they have pT > 0.0 and |eta| < 4.5 at this stage
      const Particles& neutrinos = apply<PromptFinalState>(event, "neutrinos").particlesByPt();
      FourMomentum mv;
      for (const Particle& nu : neutrinos) mv += nu.momentum();
      if (mv.pT() < 20*GeV) {
        numVetoedEvents[ind] += 1;
	numEarlyReject[ind] += 1;
	vetoEvent;
	}
      //const double cut = pT_cutoff*fcut*GeV;
 
      const Jets alljets  = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 15*GeV && Cuts::abseta < 2.5);
      Jets bjets, ljets;
      // Don't count jets that overlap with the hard leptons
      // and then tag the two bjets and remaining lighter jets
      foreach (const Jet& jet, alljets) {
        if (jet.bTagged()) {
          bjets.push_back(jet);
        } else {
          ljets.push_back(jet);
        }
      }
      //cout << "Number of b-jets = " << bjets.size() << endl;
      //cout << "Number of l-jets = " << ljets.size() << endl;
      if (bjets.size() != 2) {
        MSG_DEBUG("Event failed post-lepton-isolation b-tagging cut");
	numRejectBjets[ind] += 1;
        numVetoedEvents[ind] += 1;
        vetoEvent;
      }
      if (ljets.size() < 2) {
        MSG_DEBUG("Event failed since not enough light jets remaining after lepton-isolation");
	numRejectLjets[ind] += 1;
        numVetoedEvents[ind] += 1;
        vetoEvent;
      }
      // Boosted selection: lepton/jet overlap
      const double transmass = sqrt( 2 * lepton.pT() * mv.pT() * (1 - cos(deltaPhi(lepton, mv))) );
      if (transmass + mv.pt() <= 60*GeV) {
        numVetoedEvents[ind] += 1;
	numRejectOverlap[ind] += 1;
	vetoEvent;
      }
      // at least one b-tagged jet is required to be within R = 1 of the lepton direction
      double dist = 1.5;
      int lepJetIndex = -1;
      for (size_t i = 0; i < alljets.size(); ++i) {
        const Jet& jet = alljets[i];
        if (deltaR(jet, lepton) < dist) {
          lepJetIndex = i;
          break;
        }
      }
      if (lepJetIndex < 0){
        numVetoedEvents[ind] += 1;
	numRejectDirection[ind] += 1;
	vetoEvent;
	}
      const Jet& ljet = alljets[lepJetIndex];
	// at least one fatjet is required to be further than R = 1.5 of the lepton direction
	// and have significant angular separation 
      // Boosted selection: lepton-jet/fat-jet matching
      int fatJetIndex = -1;
      for (size_t j = 0; j < st_jets.size(); ++j) {
        const Jet& fjet = st_jets[j];
        const double dR_fatjet = deltaR(ljet, fjet);
        const double dphi = deltaPhi(lepton, fjet);
        if (dR_fatjet > 1.5 && dphi > 2.3) {
          fatJetIndex = j;
          break;
        }
      }
      if (fatJetIndex < 0){
        numVetoedEvents[ind] += 1;
	numRejectMatching[ind] += 1;
	vetoEvent;
	}
      const PseudoJet& fjet = st_jets[fatJetIndex];
      Jet jjet = fjet;

      if (ind==0) {
          numEventsPassedPtCut[ind] += 1;
	  vector<double> ratios = getRatios(fjet);
          _antikt_jet_trimmed_tauN21->fill(ratios[0],weight); 
          _antikt_jet_trimmed_tauN32->fill(ratios[1],weight);
	  _antikt_jet_trimmed_pT->fill(fjet.pt()/GeV,weight);
	  _antikt_jet_trimmed_eta->fill(fjet.eta()/GeV,weight);
	  _antikt_jet_trimmed_mass->fill(fjet.m()/GeV,weight);
      }
      if (ind==1) {
	  numEventsPassedPtCut[ind] += 1;
          vector<double> ratios = getRatios(fjet);
          _fatjet_trimmed_tauN21->fill(ratios[0],weight); 
          _fatjet_trimmed_tauN32->fill(ratios[1],weight);
          _fatjet_trimmed_eta->fill(fjet.eta()/GeV,weight);
          _fatjet_trimmed_pT->fill(fjet.pt()/GeV,weight);
          _fatjet_trimmed_mass->fill(fjet.m()/GeV,weight);
      }

      FourMomentum pbjet1; //Momentum of bjet1
      FourMomentum pbjet2; //Momentum of bjet2
      if ( deltaR(bjets[0], lepton) <= deltaR(bjets[1], lepton) ) {
        pbjet1 = bjets[0].momentum();
        pbjet2 = bjets[1].momentum();
      } else {
        pbjet1 = bjets[1].momentum();
        pbjet2 = bjets[0].momentum();
      }

      double pz = computeneutrinoz(lepton.momentum(), mv);
      FourMomentum ppseudoneutrino( sqrt(sqr(mv.px()) + sqr(mv.py()) + sqr(pz)), mv.px(), mv.py(), pz);
      FourMomentum t_leptonic = lepton.momentum() + ppseudoneutrino + pbjet1;
 
      const Jet t_lepton(t_leptonic);
      const FourMomentum t_hadronic(fjet.E(),fjet.px(),fjet.py(),fjet.pz());
      const FourMomentum Zp = t_hadronic + t_leptonic;
      const Particles& all_final_particles = apply<VetoedFinalState>(event, "VFS").particlesByPt();
      foreach( const Particle& p, all_final_particles ) {
         if (ind==0) _Et->fill(p.Et()/GeV,weight);
      }
// Now find the top and antitop at parton level
      const HepMC::GenEvent* genEvent= event.genEvent();
      for ( HepMC::GenEvent::particle_const_iterator p = genEvent->particles_begin();  p != genEvent->particles_end(); ++p ) {
        bool found_top_decay = 0;
        bool found_antitop_decay = 0;
	if((*p)->pdg_id() == 6) {	  
	  Particle top(*p);
	  HepMC::GenVertex* out_atop = (*p)->end_vertex();
	  for( HepMC::GenVertex::particle_iterator aW_finder=out_atop->particles_begin(HepMC::children) ; aW_finder != out_atop->particles_end(HepMC::children ); aW_finder++) {
      	    if ((*aW_finder)->pdg_id()== 24) {
              found_top_decay = 1;
              Tops.push_back(top);
              break;
	    }	
	  }	
	  if (found_top_decay) {
            if (ind==0) {
	      _antikt_ldiff_pT->fill((top.pT()-t_leptonic.pT())/top.pT()/GeV, weight);
	      _antikt_ldiff_eta->fill((abs(top.eta())-abs(t_leptonic.eta()))/GeV, weight);
	    }
            if (ind==1) { 
	      _fatjet_ldiff_pT->fill((top.pT()-t_leptonic.pT())/top.pT()/GeV, weight);
	      _fatjet_ldiff_eta->fill((abs(top.eta())-abs(t_leptonic.eta()))/GeV, weight);
	    }
	  }
	}  
	if((*p)->pdg_id() == -6 ) {
	  Particle top(*p);
	  HepMC::GenVertex* out_atop = (*p)->end_vertex();
	  for( HepMC::GenVertex::particle_iterator aW_finder=out_atop->particles_begin(HepMC::children) ; aW_finder != out_atop->particles_end(HepMC::children ); aW_finder++) {
            if ((*aW_finder)->pdg_id()== -24) {
              found_antitop_decay = 1;
              AntiTops.push_back(top);
	      _dR->fill(deltaR(top,jjet),weight);
              if(deltaR(top,jjet) < 1.0) {
                numCorrect[ind] += 1;
		if(ind==0) {
                  _matched_mass->fill(fjet.m()/GeV,weight);
                  _matched_diff_pT->fill((top.pT()-t_hadronic.pT())/top.pT()/GeV, weight);
                  _matched_diff_eta->fill((abs(top.eta())-abs(t_hadronic.eta()))/GeV, weight);
                  PseudoJets constits = fastjet::sorted_by_pt(fjet.constituents());
                  const double pT_ratio = constits[constits.size()-1].pt()/constits[0].pt();
                  _matched_ratio->fill(pT_ratio/GeV,weight);
                  _matched_numSJ->fill(constits.size(), weight);
                }
	        if(ind==1) {
                  _fj_matched_mass->fill(fjet.m()/GeV,weight);
                  _fj_matched_diff_pT->fill((top.pT()-t_hadronic.pT())/top.pT()/GeV, weight);
                  _fj_matched_diff_eta->fill((abs(top.eta())-abs(t_hadronic.eta()))/GeV, weight);
                  PseudoJets constits = fastjet::sorted_by_pt(fjet.constituents());
                  const double pT_ratio = constits[constits.size()-1].pt()/constits[0].pt();
                  _fj_matched_ratio->fill(pT_ratio/GeV,weight);
                  _fj_matched_numSJ->fill(constits.size(), weight);
                }
               }
               else {
                 numWrong[ind] += 1;
                }
                break;
             }
	  }
	  if (found_antitop_decay) {
            if (ind==0) {
	       _antikt_diff_pT->fill((top.pT()-t_hadronic.pT())/top.pT()/GeV, weight);
   	       _antikt_diff_eta->fill((abs(top.eta())-abs(t_hadronic.eta()))/GeV, weight);
 	    }
            if (ind==1) { 
	      _fatjet_diff_pT->fill((top.pT()-t_hadronic.pT())/top.pT()/GeV, weight);
	      _fatjet_diff_eta->fill((abs(top.eta())-abs(t_hadronic.eta()))/GeV, weight);
	    }
	    break;
	  }
        }
     }
      PseudoJets constits = fastjet::sorted_by_pt(fjet.constituents());
      const double pT_ratio = constits[constits.size()-1].pt()/constits[0].pt();
      const double sScale = constits[constits.size()-1].pt()*constits[constits.size()-1].delta_R(constits[0])/1.0; 
      if (ind==0) {
        _antikt_t_mass->fill(t_leptonic.mass()/GeV, weight);
        _antikt_t_pT->fill(t_leptonic.pT()/GeV, weight);
        _antikt_Zp_mass->fill(Zp.mass()/GeV, weight);
        _antikt_Zp_pT->fill(Zp.pT()/GeV, weight);
	_antikt_ratio->fill(pT_ratio/GeV,weight);
	_antikt_d12->fill(sScale/GeV,weight);
        _antikt_numSJ->fill(constits.size(), weight);
	}
      if (ind==1) {
        _fatjet_t_mass->fill(t_leptonic.mass()/GeV, weight);
        _fatjet_t_pT->fill(t_leptonic.pT()/GeV, weight);
        _fatjet_Zp_mass->fill(Zp.mass()/GeV, weight);
        _fatjet_Zp_pT->fill(Zp.pT()/GeV, weight);
      }
    }
/***************************************************************************************/
    // trim and recluster small jets
    void ClusterSmallJets(const Event& event, const double Rsmall, const double fcut) {
      const double weight = event.weight();
      PseudoJets small_pjets;
      PseudoJets small_ktjets;
      const Jets smalljets  = apply<FastJets>(event, "jets").jetsByPt(Cuts::pT > 15*GeV && Cuts::abseta < ymax);
      const Jets ktjets  = apply<FastJets>(event, "ktjets").jetsByPt(Cuts::abseta < ymax);
      for (const PseudoJet& pjet : smalljets) small_pjets += pjet;
      for (const PseudoJet& pjet : ktjets) small_ktjets += pjet;
// Median-area method to reduce pileup
// Jet definition needs to be the same as in the definition of ktjets
      fastjet::JetDefinition jet_def_small (fastjet::kt_algorithm, 0.5, fastjet::E_scheme, fastjet::Best);
      fastjet::Selector rho_range = fastjet::SelectorAbsRapMax(ymax);
      fastjet::JetMedianBackgroundEstimator bge(rho_range,jet_def_small,*area_def);
      bge.set_jets(small_ktjets);
      double rho = bge.rho();
      double sigma = bge.sigma();	
// Fill out subjet histograms
      _subjet_rho_distrib->fill(rho/GeV,weight);
/*
      foreach (const PseudoJet& cjet, small_ktjets) {
        _subjet_pT_distrib->fill(cjet.pt()/GeV,weight);
        _subjet_m_distrib->fill(cjet.m()/GeV,weight);
        _subjet_area_distrib->fill(cjet.area()/GeV,weight);
        if (cjet.area()!=0) _subjet_ratio_distrib->fill(cjet.pt()/cjet.area()/GeV,weight);
      }
*/
      if (rho==0) numQuietEvents += 1;
      fastjet::Subtractor subtractor(rho);
      small_pjets = subtractor(small_pjets);
// Now have antikt R = 0.3 corrected subjets which should be reclustered
// Now recluster small jets
// bottom up approach
      fastjet::JetDefinition jet_def (fastjet::antikt_algorithm, 1.0, fastjet::E_scheme, fastjet::Best);
      const fastjet::ClusterSequence kt_cs1(small_pjets, jet_def);
      PseudoJets st_jets = kt_cs1.inclusive_jets();
      PartonLevel(fastjet::sorted_by_pt(st_jets),event,0);
// Now use traditional approach
      PseudoJets fat_pjets = apply<FastJets>(event, "fat_jets").pseudoJetsByPt();
      fat_pjets = subtractor(fat_pjets);
      PseudoJets trimmed_fat_pjets;
      fastjet::Filter trimmer(fastjet::JetDefinition(fastjet::kt_algorithm, 0.3), fastjet::SelectorPtFractionMin(fcut),rho);
      for (const PseudoJet& pjet : fat_pjets) trimmed_fat_pjets += trimmer(pjet);
// Now use median method to correct the larger jets
      trimmed_fat_pjets = fastjet::sorted_by_pt(subtractor(trimmed_fat_pjets));
      PartonLevel(trimmed_fat_pjets,event,1);
   }
/***************************************************************************************/
// Unable to link Rivet with Nsubjettiness fastjet extension so just code it 
    vector<double> getRatios(PseudoJet tj) {
      double Rcut=10;
      double tau_32 = 0;
      double tau_21 = 0;
      const PseudoJets constituents = tj.constituents();
       
      if(constituents.size() < 3) tau_32=0;
      if(constituents.size() < 2) tau_21=0;
      else
        {
          PseudoJets axis3 = GetAxes(3,constituents,M_PI/2.0);
          PseudoJets axis2 = GetAxes(2,constituents,M_PI/2.0);
          PseudoJets axis1 = GetAxes(1,constituents,M_PI/2.0);
      
          double tau_3 = TauValue(1.0, constituents, axis3, Rcut) ;
          double tau_2 = TauValue(1.0, constituents, axis2, Rcut) ;
          double tau_1 = TauValue(1.0, constituents, axis1, Rcut) ;
          tau_32 =safediv( tau_3, tau_2, 0);
          tau_21 =safediv( tau_2, tau_1, 0);
        }
      vector<double> ratios;
      ratios += tau_21; 
      ratios += tau_32; 
      return ratios;
    }
/***************************************************************************************/
// Get Axes for N-subjettiness calculation
    PseudoJets GetAxes(int n_jets,const PseudoJets& constituents, double subR) const {
      if (constituents.size() < (size_t) n_jets) return constituents;
      fastjet::ClusterSequence sub_clust_seq(constituents,fastjet::JetDefinition(fastjet::kt_algorithm,subR,fastjet::E_scheme,fastjet::Best));
      return sub_clust_seq.exclusive_jets(n_jets);
    }
/***************************************************************************************/
// Calculate tau value for N-subjettiness calculation   
    double TauValue(double jet_rad,const PseudoJets& particles, const PseudoJets& axes, double Rcut) const {
      double tauNum = 0.0;
      double tauDen = 0.0;
      if (particles.size() == 0 ) return 0.0;
      for (size_t i = 0; i < particles.size(); i++) {
        double minR = 1000;
          for(size_t j = 0; j < axes.size(); j++) {
            double tempR = sqrt(particles[i].squared_distance(axes[j]));
            if (tempR < minR) minR = tempR;
          }
        if (minR > Rcut) minR = Rcut;
        tauNum += particles[i].perp() * minR;
        tauDen += particles[i].perp() * jet_rad;
      }
      return safediv(tauNum,tauDen,0);
    }
// Finished with Nsubjettiness functions
/***************************************************************************************/

    // Set up projections and book histograms
    void init() {
 
      // Base final state definition
      const FinalState fs(Cuts::abseta < 4.5);
      declare(fs, "final_state");

      // Neutrinos for MET
      IdentifiedFinalState nu_id;
      nu_id.acceptNeutrinos();
      PromptFinalState neutrinos(nu_id);
      neutrinos.acceptTauDecays(true);
      declare(neutrinos, "neutrinos");

      // Get photons used to dress leptons
      IdentifiedFinalState photons(fs);
      photons.acceptId(PID::PHOTON);

      // Use all bare muons as input to the DressedMuons projection
      IdentifiedFinalState mu_id(fs);
      mu_id.acceptIdPair(PID::MUON);
      PromptFinalState bare_mu(mu_id);
      bare_mu.acceptTauDecays(true);

      // Use all bare electrons as input to the DressedElectrons projection
      IdentifiedFinalState el_id(fs);
      el_id.acceptIdPair(PID::ELECTRON);
      PromptFinalState bare_el(el_id);
      bare_el.acceptTauDecays(true);

      // Use all bare leptons including taus for single-lepton filter
      IdentifiedFinalState lep_id(fs);
      lep_id.acceptIdPair(PID::MUON);
      lep_id.acceptIdPair(PID::ELECTRON);
      PromptFinalState bare_lep(lep_id);
      declare(bare_lep, "bare_lep");

      // Tau finding
      /// @todo Use TauFinder
      UnstableFinalState ufs;
      IdentifiedFinalState tau_id(ufs);
      tau_id.acceptIdPair(PID::TAU);
      PromptFinalState bare_tau(tau_id);
      declare(bare_tau, "bare_tau");

      // Muons and electrons must have |eta| < 2.5
      Cut eta_ranges = Cuts::abseta < 2.5;
      
      // A FinalState is used to select particles within |eta| < 2.5 and with pT
      // > 25 GeV, out of which the ChargedLeptons projection picks only the
      // electrons and muons, to be accessed later as "LFS".
      ChargedLeptons lfs(FinalState(-2.5, 2.5, 25*GeV));
      declare(lfs, "LFS");

      // Get dressed muons and the good muons (pt>25GeV)
      DressedLeptons all_dressed_mu(photons, bare_mu, 0.1, eta_ranges, true, true);
      DressedLeptons dressed_mu(photons, bare_mu, 0.1, eta_ranges & (Cuts::pT > 25.0*GeV), true, true);
      declare(dressed_mu, "muons");

      // Get dressed electrons and the good electrons (pt>25GeV)
      DressedLeptons all_dressed_el(photons, bare_el, 0.1, eta_ranges, true, true);
      DressedLeptons dressed_el(photons, bare_el, 0.1, eta_ranges & (Cuts::pT > 25.0*GeV), true, true);
      declare(dressed_el, "electrons");
      // Jet clustering

      VetoedFinalState vfs(fs);
      vfs.addVetoOnThisFinalState(lfs);
      vfs.addVetoOnThisFinalState(all_dressed_el);
      vfs.addVetoOnThisFinalState(all_dressed_mu);
      vfs.addVetoOnThisFinalState(neutrinos);
      declare(vfs, "VFS");
      // Small-R jets
      /// @todo Use extra constructor args
      //areaspec = new fastjet::GhostedAreaSpec(2.5,1,0.01);
      //area_def = new fastjet::AreaDefinition(fastjet::active_area_explicit_ghosts,*areaspec);
      FastJets jets(vfs, FastJets::ANTIKT, Rsmall);
      jets.useJetArea(area_def);
      jets.useInvisibles(JetAlg::ALL_INVISIBLES);
      jets.useMuons(JetAlg::DECAY_MUONS);
      declare(jets, "jets");

      // jets to use with median method subtraction
      /// @todo Use extra constructor args
      FastJets ktjets(vfs, FastJets::KT, 0.5);//Rsmall);
      ktjets.useJetArea(area_def);
      ktjets.useInvisibles(JetAlg::ALL_INVISIBLES);
      ktjets.useMuons(JetAlg::DECAY_MUONS);
      declare(ktjets, "ktjets");

      declare(MissingMomentum(vfs), "MissingET");

      // Large-R jets
      /// @todo Use extra constructor args
      FastJets large_jets(vfs, FastJets::ANTIKT, 1.0);
      large_jets.useJetArea(area_def);
      large_jets.useInvisibles(JetAlg::ALL_INVISIBLES);
      large_jets.useMuons(JetAlg::DECAY_MUONS);
      declare(large_jets, "fat_jets");

      // Booking of histograms
      _subjet_pT_distrib = bookHisto1D("subjet_pT",100,0.0,300);
      _subjet_m_distrib = bookHisto1D("subjet_m",40,0.0,120);
      _subjet_area_distrib = bookHisto1D("subjet_area",40,0.0,1.2);
      _subjet_rho_distrib = bookHisto1D("subjet_rho",50,0.0,200);
      _subjet_ratio_distrib = bookHisto1D("subjet_ratio",50,0.0,300);

      _fj_matched_mass = bookHisto1D("fj_matched_mass", 100, 100.0, 300.0);
      _fj_matched_diff_pT = bookHisto1D("fj_matched_diff_pT", 100, -1.0, 1.0);
      _fj_matched_diff_eta = bookHisto1D("fj_matched_diff_eta", 50, -1, 1);
      _fj_matched_ratio = bookHisto1D("fj_matched_ratio", 100, 0.0, 1.05);
      _fj_matched_numSJ = bookHisto1D("fj_matched_numSJ", 10, 0.0, 10.0);

      _matched_mass = bookHisto1D("matched_mass", 100, 100.0, 300.0);
      _matched_diff_pT = bookHisto1D("matched_diff_pT", 100, -1.0, 1.0);
      _matched_diff_eta = bookHisto1D("matched_diff_eta", 50, -1, 1);
      _matched_ratio = bookHisto1D("matched_ratio", 100, 0.0, 1.05);
      _matched_numSJ = bookHisto1D("matched_numSJ", 10, 0.0, 10.0);

      _missing_Et = bookHisto1D("missing_Et",13,0.0,520.0);
      _Et = bookHisto1D("Et",200,0.0,4000.0); 
      _dR = bookHisto1D("dR", 50, 0.0, 6.0);

      _antikt_ratio = bookHisto1D("antikt_ratio", 100, 0.0, 1.05);
      _antikt_d12 = bookHisto1D("antikt_d12", 100, 0.0, 200.0);
      _antikt_diff_pT = bookHisto1D("antikt_diff_pT", 100, -1.0, 1.0);
      _antikt_ldiff_pT = bookHisto1D("antikt_ldiff_pT", 100, -1.0, 1.0);
      _antikt_diff_eta = bookHisto1D("antikt_diff_eta", 50, -1, 1);
      _antikt_ldiff_eta = bookHisto1D("antikt_ldiff_eta", 50, -1, 1);
      _fatjet_diff_pT = bookHisto1D("fatjet_diff_pT", 100, -1.0, 1.0);
      _fatjet_ldiff_pT = bookHisto1D("fatjet_ldiff_pT", 100, -1.0, 1.0);
      _fatjet_diff_eta = bookHisto1D("fatjet_diff_eta", 50, -1, 1);
      _fatjet_ldiff_eta = bookHisto1D("fatjet_ldiff_eta", 50, -1, 1);

      //_antikt_Zp_mass = bookHisto1D("antikt_Zp_mass", 100, 500.0, 1500.0);
      _antikt_Zp_mass = bookHisto1D("antikt_Zp_mass", 200, 1000.0, 5000.0);
      _antikt_Zp_pT = bookHisto1D("antikt_Zp_pT", 30, 0.0, 300.0);
      _antikt_t_mass = bookHisto1D("antikt_t_leptonic_mass", 100, 100.0, 300.0);
      _antikt_t_pT = bookHisto1D("antikt_t_leptonic_pT", 30, 0.0, 600.0);
      _antikt_jet_trimmed_mass = bookHisto1D("antikt_jet_trimmed_mass", 100, 100.0, 300.0);
      _antikt_jet_trimmed_pT = bookHisto1D("antikt_jet_trimmed_pT", 100, 350.0, 2000.0);
      _antikt_jet_trimmed_eta = bookHisto1D("antikt_jet_trimmed_eta",20 ,0,2.5);
      _antikt_jet_trimmed_tauN32 = bookHisto1D("antikt_jet_trimmed_tauN32", 20, 0.0, 1.2);
      _antikt_jet_trimmed_tauN21 = bookHisto1D("antikt_jet_trimmed_tauN21", 20, 0.0, 1.2);
      _antikt_numSJ = bookHisto1D("numSJ", 10, 0.0, 10.0);

      //_fatjet_Zp_mass = bookHisto1D("fatjet_Zp_mass", 100, 500.0, 1500.0);
      _fatjet_Zp_mass = bookHisto1D("fatjet_Zp_mass", 200, 1000.0, 5000.0);
      _fatjet_Zp_pT = bookHisto1D("fatjet_Zp_pT", 30, 0.0, 300.0);
      _fatjet_t_mass = bookHisto1D("fatjet_t_leptonic_mass", 100, 100.0, 300.0);
      _fatjet_t_pT = bookHisto1D("fatjet_t_leptonic_pT", 30, 0.0, 600.0);
      _fatjet_trimmed_pT = bookHisto1D("fatjet_trimmed_pT", 100, 350.0, 2000.0);
      _fatjet_trimmed_eta = bookHisto1D("fatjet_trimmed_eta",20 ,0,2.5);
      _fatjet_trimmed_mass = bookHisto1D("fatjet_trimmed_mass", 100, 100.0, 300.0);
      _fatjet_trimmed_d12 = bookHisto1D("fatjet_trimmed_d12", 20, 0.0, 200.0);
      _fatjet_trimmed_tauN21 = bookHisto1D("fatjet_trimmed_tauN21", 20, 0.0, 1.2);
      _fatjet_trimmed_tauN32 = bookHisto1D("fatjet_trimmed_tauN32", 20, 0.0, 1.2);
   }
/***************************************************************************************/

// Begin analysis
    void analyze(const Event& event) {
      const double weight = event.weight();
      numEvents += 1;
      const MissingMomentum& met = apply<MissingMomentum>(event, "MissingET");
      _missing_Et->fill(met.vectorEt().mod()/GeV,weight);
      ClusterSmallJets(event,Rsmall,fcut); 
    }
/***************************************************************************************/
    void finalize() {
     for (int i=0;i<2;i++) {
        cout << "Number of Events = " << numEvents << endl;
        cout << "Number of Quiet Events = " << numQuietEvents << endl;
        cout << "Number of Early rejections from leptons and neutrinos = " << numEarlyReject[i] << endl;
        cout << "Number of rejections for not enough b jets = " << numRejectBjets[i] << endl;
        cout << "Number of rejections for not enough l jets = " << numRejectLjets[i] << endl;
        cout << "Number of rejections for no overlap = " << numRejectOverlap[i] << endl;
        cout << "Number of rejections for direction = " << numRejectDirection[i] << endl;
        cout << "Number of rejections for no matching = " << numRejectMatching[i] << endl;
        cout << "Number of Vetoed Events = " << numVetoedEvents[i] << endl;
        cout << "Number of Events with Pt > 350 GeV = " << numEventsPassedPtCut[i] << endl;
        if (numCorrect[i] != 0 || numWrong[i] != 0) cout << "Percent of Jets chosen (for method " << i << " ) that correctly have a top quark inside = " << numCorrect[i]/(numCorrect[i]+numWrong[i]) * 100 << " %" << endl;
        cout << "Number of Correctly Identified tops: " << numCorrect[i] << endl;
        cout << "Number of Incorrectly Identified tops: " << numWrong[i] << endl;
      }
      cout << "Scale Factor = " << 2*crossSection()*L/femtobarn / sumOfWeights() << endl;
      cout << "Cross Section = " << crossSection() << endl;
      cout << "femtobarn = " << femtobarn << endl;
      cout << "sumOfWeights = " << sumOfWeights() << endl;
      scale(_subjet_pT_distrib, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_subjet_rho_distrib, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_subjet_m_distrib, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_subjet_area_distrib, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_subjet_ratio_distrib, 2*crossSection()*L/femtobarn / sumOfWeights());

      scale(_fj_matched_ratio, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fj_matched_diff_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fj_matched_numSJ, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fj_matched_diff_eta, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fj_matched_mass, 2*crossSection()*L/femtobarn / sumOfWeights());

      scale(_matched_ratio, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_matched_diff_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_matched_numSJ, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_matched_diff_eta, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_matched_mass, 2*crossSection()*L/femtobarn / sumOfWeights());
     
      scale(_dR, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_ratio, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_d12, 2*crossSection()*L/femtobarn / sumOfWeights());
 
      scale(_antikt_diff_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_ldiff_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_diff_eta, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_ldiff_eta, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_diff_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_ldiff_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_diff_eta, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_ldiff_eta, 2*crossSection()*L/femtobarn / sumOfWeights());
     
      scale(_missing_Et, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_Et, 2*crossSection()*L/femtobarn / sumOfWeights());
      //scale(_missing_Et, 2*crossSection()*L/femtobarn / sumOfWeights());

      scale(_antikt_Zp_mass, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_Zp_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_t_mass, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_t_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_jet_trimmed_mass, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_jet_trimmed_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_jet_trimmed_eta, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_jet_trimmed_tauN32, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_antikt_jet_trimmed_tauN21, 2*crossSection()*L/femtobarn / sumOfWeights());

      scale(_fatjet_Zp_mass, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_Zp_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_t_mass, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_t_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_trimmed_mass, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_trimmed_pT, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_trimmed_eta, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_trimmed_d12, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_trimmed_tauN21, 2*crossSection()*L/femtobarn / sumOfWeights());
      scale(_fatjet_trimmed_tauN32, 2*crossSection()*L/femtobarn / sumOfWeights());
    }

  private:

    double ymax= 2.5; // e.g. if particles go up to y=2.5
    int nrepeat = 1;
    double ghost_area = 0.01;
    fastjet::GhostedAreaSpec* areaspec;//(ymax, nrepeat, ghost_area);
    fastjet::AreaDefinition* area_def;//(fastjet::active_area_explicit_ghosts, areaspec);
 
    Particles Tops;
    Particles AntiTops;
// Bunch of things to keep track of how many vetoes from each vetoEvent
    int numEvents = 0;
    int numQuietEvents = 0;
    int numEmpty[2] = {0,0};
    int numEarlyReject[2] = {0,0};
    int numRejectBjets[2] = {0,0};
    int numRejectLjets[2] = {0,0};
    int numVetoedEvents[2] = {0,0};
    double numCorrect[2] = {0,0};
    double numWrong[2] = {0,0};
    int numEventsPassedPtCut[2] = {0,0};
    int numRejectOverlap[2] = {0,0};
    int numRejectDirection[2] = {0,0};
    int numRejectMatching[2] = {0,0};
// Integrated Luminosity
    double L = 20.3;
// Parameters for small-R trimming
    const double Rsmall = 0.3, fcut = 0.1;
// pT cutoff for R=1.0 jets (reclustered and large-R)
    const double pT_cutoff = 350;
    const double mass_cutoff = 100;

    Histo1DPtr _fj_matched_mass;
    Histo1DPtr _fj_matched_diff_pT;
    Histo1DPtr _fj_matched_diff_eta;
    Histo1DPtr _fj_matched_ratio;
    Histo1DPtr _fj_matched_numSJ;

    Histo1DPtr _matched_mass;
    Histo1DPtr _matched_diff_pT;
    Histo1DPtr _matched_diff_eta;
    Histo1DPtr _matched_ratio;
    Histo1DPtr _matched_numSJ;

    Histo1DPtr _dR;
    Histo1DPtr _antikt_ratio;
    Histo1DPtr _antikt_d12;
    Histo1DPtr _antikt_diff_pT;
    Histo1DPtr _antikt_ldiff_pT;
    Histo1DPtr _antikt_diff_eta;
    Histo1DPtr _antikt_ldiff_eta;
    Histo1DPtr _fatjet_diff_pT;
    Histo1DPtr _fatjet_ldiff_pT;
    Histo1DPtr _fatjet_diff_eta;
    Histo1DPtr _fatjet_ldiff_eta;

    Histo1DPtr _antikt_Zp_mass;
    Histo1DPtr _antikt_Zp_pT;
    Histo1DPtr _antikt_t_mass;
    Histo1DPtr _antikt_t_pT;
    Histo1DPtr _antikt_jet_trimmed_mass;
    Histo1DPtr _antikt_jet_trimmed_pT;
    Histo1DPtr _antikt_jet_trimmed_eta;
    Histo1DPtr _antikt_jet_trimmed_tauN32;
    Histo1DPtr _antikt_jet_trimmed_tauN21;

    Histo1DPtr _fatjet_Zp_mass;
    Histo1DPtr _fatjet_Zp_pT;
    Histo1DPtr _fatjet_t_mass;
    Histo1DPtr _fatjet_t_pT;
    Histo1DPtr _fatjet_trimmed_mass;
    Histo1DPtr _fatjet_trimmed_pT;
    Histo1DPtr _fatjet_trimmed_eta;
    Histo1DPtr _fatjet_trimmed_d12;
    Histo1DPtr _fatjet_trimmed_tauN21;
    Histo1DPtr _fatjet_trimmed_tauN32;

    Histo1DPtr _missing_Et;
    Histo1DPtr _Et;
    Histo1DPtr _antikt_numSJ;
    Histo1DPtr _subjet_pT_distrib;
    Histo1DPtr _subjet_rho_distrib;
    Histo1DPtr _subjet_m_distrib;
    Histo1DPtr _subjet_area_distrib;
    Histo1DPtr _subjet_ratio_distrib;
  };

  // The hook for the plugin system
  DECLARE_RIVET_PLUGIN(tt_signal);
}