// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Projections/DressedLeptons.hh"
#include "Rivet/Projections/LeadingParticlesFinalState.hh"
#include "Rivet/Projections/MissingMomentum.hh"
#include "Rivet/Projections/WFinder.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Projections/FinalState.hh"

namespace Rivet {


  class WFINDER_TEST : public Analysis {
  public:

    /// Constructor
    WFINDER_TEST()
      : Analysis("WFINDER_TEST")
    {
      _dR = 0.2;
      _lepton=PID::MUON;
    }


    /// Book histograms and initialise projections before the run
    void init() {

      FinalState fs;
      IdentifiedFinalState allleptons;
      // allleptons.acceptIdPair(PID::ELECTRON);
      allleptons.acceptIdPair(PID::MUON);
      Cut cuts = Cuts::abseta < 2.5 && Cuts::pT > 20*GeV;
      DressedLeptons leptons(fs, allleptons, 0.1, cuts);
      addProjection(leptons, "leptons");

      // Leading neutrinos for Etmiss
      LeadingParticlesFinalState neutrinos(fs);
      // neutrinos.addParticleIdPair(PID::NU_E);
      neutrinos.addParticleIdPair(PID::NU_MU);
      neutrinos.setLeadingOnly(true);
      addProjection(neutrinos, "neutrinos");
      
      WFinder wfinder(fs, Cuts::abseta < 3.0 && Cuts::pT > 20*GeV, _lepton, 0.0*GeV, 1000.0*GeV, 25.0*GeV, _dR);
      addProjection(wfinder, "WFinder");

      _h_W_mass = bookHisto1D("W_mass", 100, 0, 300);
      _h_W_massCheck = bookHisto1D("W_massCheck", 100, 0, 300);

      
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {
      const double weight = event.weight();

      const WFinder& wfinder = applyProjection<WFinder>(event, "WFinder");
      if (wfinder.bosons().size() != 1) {
        vetoEvent;
      }
      FourMomentum wmom(wfinder.bosons().front().momentum());
      _h_W_mass->fill(wmom.mass()/GeV, weight);

      // implementation of simple Wfinder 
      const vector<DressedLepton>& leptons = applyProjection<DressedLeptons>(event, "leptons").dressedLeptons();
      Particles neutrinos = applyProjection<FinalState>(event, "neutrinos").particlesByPt();
      if (leptons.size() != 1 || (neutrinos.size() == 0)) {
	vetoEvent;
      }

      FourMomentum lepton = leptons[0].momentum();
      FourMomentum p_miss = neutrinos[0].momentum();
      if (p_miss.Et() < 25.0*GeV) {
	vetoEvent;
      }

      double invMass = sqrt(2.0 * lepton.pT()/GeV * p_miss.pT()/GeV * (cosh(lepton.eta()-p_miss.eta()) - cos( lepton.phi()-p_miss.phi()) ) );
      _h_W_massCheck->fill(invMass, weight);
      
    }


    /// Normalise histograms etc., after the run
    void finalize() {

      scale(_h_W_mass, crossSection()/picobarn/sumOfWeights());
      scale(_h_W_massCheck, crossSection()/picobarn/sumOfWeights());

    }


  private:

    double _dR;
    PdgId _lepton;

    Histo1DPtr _h_W_mass;
    Histo1DPtr _h_W_massCheck;
    
  };



  // The hook for the plugin system
  DECLARE_RIVET_PLUGIN(WFINDER_TEST);


}