// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/RivetAIDA.hh"
#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"

namespace Rivet {


  class CMS_2011_I954992 : public Analysis {
  public:

    /// Constructor
    CMS_2011_I954992()
      : Analysis("CMS_2011_I954992"), 
	a(0.),
	b(0.),
	c(0.),
	d(0.),
	e(0.),
	f(0.),
	g(0.)
    {    }


  public:

    void init() {

      /// Get muons which pass the initial kinematic cuts
      IdentifiedFinalState muon_fs(-2.1, 2.1, 4.0*GeV);
      //IdentifiedFinalState muon_fs(-1000., 1000., 4.0*GeV);
      muon_fs.acceptIdPair(MUON);
      addProjection(muon_fs, "MUON_FS");

      _h_sigma = bookHistogram1D(1,1,1);
      h1 = bookHistogram1D("h1", 40, 0., 40., "Single muon pT, after all selections", "$p_{T}\\mu^{+}$ [GeV]", "Events/1 GeV"); // pt antimuon
      h2 = bookHistogram1D("h2", 40, 0., 40., "Single muon pT, after all selections", "$p_{T}\\mu^{-}$ [GeV]", "Events/1 GeV"); // pt muon
      h3 = bookHistogram1D("h3", 16, -2.4, 2.4, "Single muon pseudorapidity, after all selections", "$\\eta \\mu^{+}$ [GeV]", "Events/0.3 GeV"); // eta antimuon 
      h4 = bookHistogram1D("h4", 16, -2.4, 2.4, "Single muon pseudorapidity, after all selections", "$\\eta \\mu^{+}$ [GeV]", "Events/0.3 GeV"); // eta muon
      h5 = bookHistogram1D("h5", 12, -0.1, 1.1, "Delta pT(mumu), after all selections", "$fabs(\\Delta p_{T}(\\mu\\mu))$", "Events/0.1 GeV"); //deltapt
      h6 = bookHistogram1D("h6", 100, 0., 150., "dimuon mass, after all selections", "$\\mu\\mu mass$ [GeV]", "Events/1.5 GeV"); // dimuon mass
      h7 = bookHistogram1D("h7", 22, 0., 0.11, "phi, after all selections", "$\\mu\\mu 1-fabs(\\Delta\\phi / pi)$", "Events/0.005 GeV"); // dphi/pi

    }

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

      const ParticleVector& muonFS = applyProjection<IdentifiedFinalState>(event, "MUON_FS").particles();

      MSG_DEBUG("size = " << muonFS.size());

      a += weight;

      if(muonFS.size() != 2) vetoEvent;

      b += weight;

      if(PID::charge(muonFS.at(0)) != PID::charge(muonFS.at(1))) {

      c += weight;

         const double dimuon_mass = (muonFS.at(0).momentum() + muonFS.at(1).momentum()).mass();
	 const double v_angle = muonFS.at(0).momentum().angle(muonFS.at(1).momentum());
         const double dPhi = deltaPhi(muonFS.at(0).momentum().phi(), muonFS.at(1).momentum().phi());
         const double deltaPt = fabs(muonFS.at(0).momentum().pT() - muonFS.at(1).momentum().pT());

	 MSG_DEBUG("dimuon_mass = " << dimuon_mass << " dPhi = " << dPhi << " (1-fabs(dPhi/PI)) = " << (1-fabs(dPhi/PI)) << " pT1 = " << muonFS.at(0).momentum().pT()   
				    << " pT2 = " << muonFS.at(1).momentum().pT() << " deltapt = " << deltaPt );
	 
         if (dimuon_mass >= 11.5*GeV) { 
	 d += weight;
	    if (v_angle < 0.95*PI) { 
	    e += weight;
	       if ( (1-fabs(dPhi/PI)) < 0.1) {
	       f += weight;
	          if (deltaPt < 1.) {
	          g += weight;
	             _h_sigma->fill(7000/GeV, weight);
		     h5->fill(deltaPt, weight);
		     h6->fill(dimuon_mass, weight);
		     h7->fill((1-fabs(dPhi/PI)), weight);

		     if (PID::charge(muonFS.at(0))>0) { // antimuon
			h1->fill(muonFS.at(0).momentum().pT(), weight);
			h3->fill(muonFS.at(0).momentum().eta(), weight);
		     } else { // muon
			h2->fill(muonFS.at(0).momentum().pT(), weight);
			h4->fill(muonFS.at(0).momentum().eta(), weight);
		     }

		     if (PID::charge(muonFS.at(1))>0) { // antimuon
			h1->fill(muonFS.at(1).momentum().pT(), weight);
			h3->fill(muonFS.at(1).momentum().eta(), weight);
		     } else { // muon
			h2->fill(muonFS.at(1).momentum().pT(), weight);
			h4->fill(muonFS.at(1).momentum().eta(), weight);
		     }

		  }
	       }
	    }
         }
      }

    }

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

      scale(_h_sigma, crossSection()/picobarn/sumOfWeights());
      MSG_INFO("sumOfWeights, a(all), b(size), c(charge), d(mass), e(3D angle), f(phi), g(deltapt) = " << sumOfWeights() << ", " << a << ", " << b << ", " << c << ", " << d << ", " << e << ", " << f << ", " << g);
/*
      scale(h1, 1./g);
      scale(h2, 1./g);
      scale(h3, 1./g);
      scale(h4, 1./g);
      scale(h5, 1./g);
      scale(h6, 1./g);
      scale(h7, 1./g);
*/
    }

  private:

    AIDA::IHistogram1D * _h_sigma;
    AIDA::IHistogram1D * h1;
    AIDA::IHistogram1D * h2;
    AIDA::IHistogram1D * h3;
    AIDA::IHistogram1D * h4;
    AIDA::IHistogram1D * h5;
    AIDA::IHistogram1D * h6;
    AIDA::IHistogram1D * h7;

    double a;
    double b;
    double c;
    double d;
    double e;
    double f;
    double g;

  };


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

}