// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/RivetAIDA.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/Projections/UnstableFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"

namespace Rivet {

  class CMS_QCD_11_010 : public Analysis {
  public:

    CMS_QCD_11_010()
      : Analysis("CMS_QCD_11_010")
    {  }

    void init() {

      ChargedFinalState bscm(-4.65,-3.23, 0.0*GeV);
      ChargedFinalState bscp( 3.23, 4.65, 0.0*GeV);
      addProjection(bscm, "BSCM");
      addProjection(bscp, "BSCP");

      const ChargedFinalState cfs(-2.5, 2.5, 0.5*GeV); // tracks are accepted with pT > 0.5 GeV
      addProjection(cfs, "CFS");

      const FastJets jetpro(cfs, FastJets::ANTIKT, 0.5);
      addProjection(jetpro, "Jets");

      const UnstableFinalState ufs(-2.0, 2.0, 0.5*GeV);
      addProjection(ufs, "UFS");

      _hist_nTrans_Lambda = bookProfile1D(1, 1, 1);
      _hist_nTrans_Kaon = bookProfile1D(2, 1, 1);
      _hist_ptsumTrans_Lambda = bookProfile1D(3, 1, 1);
      _hist_ptsumTrans_Kaon = bookProfile1D(4, 1, 1);

    }

    void analyze(const Event& event) {

      const double weight = event.weight();

      const ChargedFinalState& bscm = applyProjection<ChargedFinalState>(event, "BSCM");
      const ChargedFinalState& bscp = applyProjection<ChargedFinalState>(event, "BSCP");
      if (bscm.size()<1 || bscp.size()<1) {
        vetoEvent;
      }

      Jets jets = applyProjection<FastJets>(event, "Jets").jetsByPt();
      if (jets.size() < 1) {
	vetoEvent;
      }
      if (fabs(jets[0].momentum().eta()) >= 2 || jets[0].momentum().pT() < 1.) {
	vetoEvent;
      }

      FourMomentum p_lead = jets[0].momentum();
      const double philead = p_lead.phi();
      const double pTlead  = p_lead.perp();

      const UnstableFinalState& ufs = applyProjection<UnstableFinalState>(event, "UFS");

      double numTrans_Kaon06(0.0), ptSumTrans_Kaon06(0.0);
      double numTrans_Lambda(0.), ptSumTrans_Lambda(0.0);

      foreach (const Particle& p, ufs.particles()) {
        double dphi = fabs(deltaPhi(philead, p.momentum().phi()));
        double eta = fabs(p.momentum().eta());
        double pT = p.momentum().pT();
        const PdgId id = p.pdgId();
      
        switch (id) {
	case 310:
	if (((PI/3.0 < dphi) && (dphi < 2.0*PI/3.0)))   //&& (pT > 0.6))
	{
	      if (pT > 0.6)
	      {
			ptSumTrans_Kaon06 += pT;	
			numTrans_Kaon06 += 1.0;
	      } 
	}
	break;

	case 3122: case -3122:
	if (((PI/3.0 < dphi) && (dphi < 2.0*PI/3.0)) && (pT > 1.5)) 
	{
              ptSumTrans_Lambda += pT;	
              numTrans_Lambda += 1.0;
	}
	break;
	}

      }

      _hist_nTrans_Kaon->fill(pTlead/GeV, numTrans_Kaon06 / ((8.0 * PI/3.0)), weight);
      _hist_nTrans_Lambda->fill(pTlead/GeV, numTrans_Lambda / ((8.0 * PI/3.0)), weight);
      _hist_ptsumTrans_Kaon->fill(pTlead/GeV , ptSumTrans_Kaon06 / ((GeV * (8.0 * PI/3.0))) , weight);
      _hist_ptsumTrans_Lambda->fill(pTlead/GeV , ptSumTrans_Lambda / ((GeV * (8.0 * PI/3.0))) , weight);

    } // end of events

    void finalize() {

    }

  private:

    AIDA::IProfile1D *_hist_nTrans_Kaon;
    AIDA::IProfile1D *_hist_nTrans_Lambda;
    AIDA::IProfile1D *_hist_ptsumTrans_Kaon;
    AIDA::IProfile1D *_hist_ptsumTrans_Lambda;

  };

  DECLARE_RIVET_PLUGIN(CMS_QCD_11_010);

}