// -*- C++ -*-
#include <iostream>
#include <sstream>
#include <string>

#include "Rivet/Analysis.hh"
#include "Rivet/RivetAIDA.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/FinalState.hh"

#include "Rivet/Projections/IdentifiedFinalState.hh"
#include "Rivet/Jet.hh"
#include "Rivet/Projections/FastJets.hh"

#include "fastjet/internal/base.hh"
#include "fastjet/JetDefinition.hh"
#include "fastjet/AreaDefinition.hh"
#include "fastjet/ClusterSequence.hh"
#include "fastjet/ClusterSequenceArea.hh"
#include "fastjet/PseudoJet.hh"


namespace Rivet {

  /// @brief Measurement of isolated diphoton + X differential cross-sections
  ///
  /// Inclusive isolated gamma gamma cross-sections, differential in M(gg), pT(gg), 
  /// dphi(gg), cos(theta*)_CS
  ///
  /// @author Giovanni Marchiori

  class ATLAS_2012_I1199269 : public Analysis {
  public:

    /// Constructor
    ATLAS_2012_I1199269()
      : Analysis("ATLAS_2012_I1199269")
    {
      setBeams(PROTON, PROTON);
      setNeedsCrossSection(true);

      _eta_bins_areaoffset.push_back(0.0);
      _eta_bins_areaoffset.push_back(1.5);
      _eta_bins_areaoffset.push_back(3.0);
    }

  public:

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

      FinalState fs;
      addProjection(fs, "FS");

      FastJets fj(fs, FastJets::KT, 0.5);
      _area_def = new fastjet::AreaDefinition(fastjet::VoronoiAreaSpec());
      fj.useJetArea(_area_def);
      addProjection(fj, "KtJetsD05");

      IdentifiedFinalState photonfs(-2.37, 2.37, 22.0*GeV);
      photonfs.acceptId(PHOTON);
      addProjection(photonfs, "Photon");

      _h_M            = bookHistogram1D(1, 1, 1); 
      _h_pT           = bookHistogram1D(2, 1, 1);
      _h_dPhi         = bookHistogram1D(3, 1, 1);
      _h_cosThetaStar = bookHistogram1D(4, 1, 1);
    }


    int getEtaBin(double eta_w, bool area_eta) const {

      double eta = fabs(eta_w);

      int v_iter = 0;

      if (!area_eta) {
	for (v_iter=0; v_iter < (int)_eta_bins.size()-1; v_iter++) {
	  if (eta >= _eta_bins.at(v_iter) && eta < _eta_bins.at(v_iter+1))
	    break;
	}
      }
      else {
	for (v_iter=0; v_iter < (int)_eta_bins_areaoffset.size()-1; v_iter++) {
	  if (eta >= _eta_bins_areaoffset.at(v_iter) && eta < _eta_bins_areaoffset.at(v_iter+1))
	    break;
	}
      }
      
      return v_iter;
    }


    /// Perform the per-event analysis
    void analyze(const Event& event) {

      const double weight = event.weight();
      
      //
      // require at least 2 photons in final state
      //
      ParticleVector photons = applyProjection<IdentifiedFinalState>(event, "Photon").particlesByPt();
      if (photons.size() < 2) {
        vetoEvent;
      }

      //
      // compute the median energy density
      //
      _ptDensity.clear();
      _sigma.clear();
      _Njets.clear();
      std::vector< std::vector<double> > ptDensities;
      std::vector<double> emptyVec;
      ptDensities.assign(_eta_bins_areaoffset.size()-1,emptyVec);
      
      const fastjet::ClusterSequenceArea* clust_seq_area = applyProjection<FastJets>(event, "KtJetsD05").clusterSeqArea();
      foreach (const fastjet::PseudoJet& jet, applyProjection<FastJets>(event, "KtJetsD05").pseudoJets(0.0*GeV)) {
	double eta = fabs(jet.eta());
	double pt = fabs(jet.perp());
	
	// get the cluster sequence
	double area = clust_seq_area->area(jet);
	
	if (area > 10e-4 && fabs(eta)<_eta_bins_areaoffset[_eta_bins_areaoffset.size()-1]) {
	  ptDensities.at(getEtaBin(fabs(eta),true)).push_back(pt/area);
	}
      }
      
      for (int b=0; b<(int)_eta_bins_areaoffset.size()-1;b++) {
	double median = 0.0;
	double sigma = 0.0;
	int Njets = 0;

	if (ptDensities[b].size() > 0) {

	  std::sort(ptDensities[b].begin(), ptDensities[b].end());
	  int nDens = ptDensities[b].size();
	  if( nDens % 2 == 0 )
	    median = (ptDensities[b][nDens/2]+ptDensities[b][(nDens-2)/2])/2;
	  else
	    median = ptDensities[b][(nDens-1)/2];
	  sigma = ptDensities[b][(int)(.15865*nDens)];
	  Njets = nDens;
	}
	_ptDensity.push_back(median);
	_sigma.push_back(sigma);
	_Njets.push_back(Njets);
      }


      //
      // Loop over photons and fill vector of isolated ones
      //
      ParticleVector isolated_photons;
      foreach (const Particle& photon, photons) {

	//
	// remove photons in crack
	//
        double eta_P = photon.momentum().eta();
	if (fabs(eta_P)>=1.37 && fabs(eta_P)<1.52) continue;
	
        double phi_P = photon.momentum().phi();
	
	//
	// compute isolation
	//
	
	// std EtCone
	ParticleVector fs = applyProjection<FinalState>(event, "FS").particles();
	FourMomentum mom_in_EtCone;
	foreach (const Particle& p, fs) {
	  // check if it's in the cone of .4
	  if (deltaR(eta_P, phi_P, p.momentum().eta(), p.momentum().phi()) >= 0.4) continue;
	  
	  // check if it's in the 5x7 central core
	  if (fabs(eta_P-p.momentum().eta()) < .025*7.0*0.5 && 
	      fabs(phi_P-p.momentum().phi()) < (PI/128.)*5.0*0.5) continue;
	  
	  mom_in_EtCone += p.momentum();
	}
	
	// now figure out the correction (area*density)
	float EtCone_area = PI*.4*.4 - (7.0*.025)*(5.0*PI/128.);
	float correction = _ptDensity[getEtaBin(eta_P,true)]*EtCone_area;
	
      	// shouldn't need to subtract photon
	// note: using expected cut at hadron/particle level, not cut at reco level
	if (mom_in_EtCone.Et()-correction > 4.0*GeV) { 
	  continue;
	}
	
	// add photon to list of isolated ones
	isolated_photons.push_back(photon);
      }

      //
      // require at least two isolated photons
      //
      if (isolated_photons.size() < 2) {
        vetoEvent;
      }

      //
      // select leading pT pair 
      //
      std::sort(isolated_photons.begin(), isolated_photons.end(), cmpParticleByPt);
      FourMomentum y1 = isolated_photons[0].momentum();
      FourMomentum y2 = isolated_photons[1].momentum();

      //
      // leading photon should have pT>25 GeV
      //
      if (y1.pT()/GeV < 25) {
        vetoEvent;
      }
      
      //
      // require the two photons to be separated (dR>0.4)
      //
      if (deltaR(y1, y2)<0.4) {
        vetoEvent;
      }

      FourMomentum yy=y1+y2;
      double Myy = yy.mass()/GeV;
      /*
      if (Myy<30.0 || Myy>350.0) {
        vetoEvent;
      }
      */

      double pTyy = yy.pT()/GeV;
      /*
      if (Myy<pTyy) {
        vetoEvent;
      }
      */

      double dPhiyy = mapAngle0ToPi(y1.phi()-y2.phi());
      /*
      if (dPhiyy<0.5*M_PI) {
        vetoEvent;
      }
      */

      /*
	double Yyy = atanh( (y1.pT()*sinh(y1.eta())+y2.pT()*sinh(y2.eta())) / 
	(y1.pT()*cosh(y1.eta())+y2.pT()*cosh(y2.eta())) );
      */

      double costhetayy = 2*(y1.pT()/GeV)*(y2.pT()/GeV)*sinh(y1.eta()-y2.eta())/(Myy*sqrt(Myy*Myy+pTyy*pTyy));
      
      _h_M->fill(Myy, weight);
      _h_pT->fill(pTyy, weight);
      _h_dPhi->fill(dPhiyy, weight);
      _h_cosThetaStar->fill(costhetayy, weight);

    }


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

      scale(_h_M, crossSection()/sumOfWeights());
      scale(_h_pT, crossSection()/sumOfWeights());
      scale(_h_dPhi, crossSection()/sumOfWeights());
      scale(_h_cosThetaStar, crossSection()/sumOfWeights());
    }

  private:

    AIDA::IHistogram1D *_h_M;
    AIDA::IHistogram1D *_h_pT;
    AIDA::IHistogram1D *_h_dPhi;
    AIDA::IHistogram1D *_h_cosThetaStar;

    fastjet::AreaDefinition* _area_def;

    std::vector<float> _eta_bins;
    std::vector<float> _eta_bins_areaoffset;

    std::vector<float> _ptDensity;
    std::vector<float> _sigma;
    std::vector<float> _Njets;
  };


  DECLARE_RIVET_PLUGIN(ATLAS_2012_I1199269);
}