// -*- C++ -*-
#include "Rivet/Analysis.hh"
#include "Rivet/RivetAIDA.hh"
#include "Rivet/Tools/Logging.hh"
#include "Rivet/Tools/BinnedHistogram.hh"
#include "Rivet/RivetBoost.hh"
#include "Rivet/Projections/DISKinematics.hh"
#include "FinalStateBreit.hh"
#include "Rivet/AnalysisLoader.hh"
#include "Rivet/Math/Constants.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
#include "InvariantJetMass.hh"
#include "Cparam.hh"
#include "Rivet/Projections/Thrust.hh"
#include "Broadening.hh"
#include "ThrustGamma.hh"
#include "BroadeningGamma.hh"

namespace Rivet {
	
  ///@brief ZEUS Event shapes in deep inelastic scattering at HERA
  ///@author Alexander Radovic edit by Jack Pegler
  ///based on ZEUS collaboration paper DESY-06-042

  class ZEUS_2007_S6591329 : public Analysis {

  public:

    /// Constructor
    ZEUS_2007_S6591329() : Analysis("ZEUS_2007_S6591329"){

      setBeams(POSITRON,PROTON);
    }



  public:

    void init() {

      /// Get Projections:
      const FinalState fs;
      addProjection(fs,"FS");

      const DISKinematics dk;
      addProjection(dk, "DK");

      const FinalStateBreit fsb= FinalStateBreit(dk,1,.25);
      addProjection(fsb, "FSB");

      const InvariantJetMass invariantjetmassb(fsb);
      addProjection(invariantjetmassb, "IJM");

      const Cparam cparamb(fsb);
      addProjection(cparamb, "CPM");

      const Thrust thrustb(fsb);
      addProjection(thrustb, "Thrust");

      const Broadening broadeningb(fsb);
      addProjection(broadeningb, "Broadening");

      const ThrustGamma thrustgammab(dk, fsb);
      addProjection(thrustgammab, "ThrustGamma");

      const BroadeningGamma broadeninggammab(dk, fsb);
      addProjection(broadeninggammab, "BroadeningGamma");


      ///book histograms:
      ///DIS Variables
      _histQ = bookHistogram1D("Q2",100, 0, 20480);
      _histX = bookHistogram1D("X",100, 0,1);
      _histY = bookHistogram1D("Y",100, 0,1);

      ///Invariant Jet Mass
      _difM.addHistogram(80.0, 160.0, bookHistogram1D(7,1,1));
      _difM.addHistogram(160.0, 320.0, bookHistogram1D(8,1,1));
      _difM.addHistogram(320.0, 640.0, bookHistogram1D(9,1,1));
      _difM.addHistogram(640.0, 1280.0, bookHistogram1D(10,1,1));
      _difM.addHistogram(1280.0, 2560.0, bookHistogram1D(11,1,1));
      _difM.addHistogram(2560.0, 5120.0, bookHistogram1D(12,1,1));
      _difM.addHistogram(5120.0, 10240.0, bookHistogram1D(13,1,1));
      _difM.addHistogram(10240.0, 20480.0, bookHistogram1D(14,1,1));
 

      ///C-Paramater
      _difC.addHistogram(80.0, 160.0, bookHistogram1D(15,1,1));
      _difC.addHistogram(160.0, 320.0, bookHistogram1D(16,1,1));
      _difC.addHistogram(320.0, 640.0, bookHistogram1D(17,1,1));
      _difC.addHistogram(640.0, 1280.0, bookHistogram1D(18,1,1));
      _difC.addHistogram(1280.0, 2560.0, bookHistogram1D(19,1,1));
      _difC.addHistogram(2560.0, 5120.0, bookHistogram1D(20,1,1));
      _difC.addHistogram(5120.0, 10240.0, bookHistogram1D(21,1,1));
      _difC.addHistogram(10240.0, 20480.0, bookHistogram1D(22,1,1));
 

      ///Thrust
      _difT.addHistogram(80.0, 160.0, bookHistogram1D(23,1,1));
      _difT.addHistogram(160.0, 320.0, bookHistogram1D(24,1,1));
      _difT.addHistogram(320.0, 640.0, bookHistogram1D(25,1,1));
      _difT.addHistogram(640.0, 1280.0, bookHistogram1D(26,1,1));
      _difT.addHistogram(1280.0, 2560.0, bookHistogram1D(27,1,1));
      _difT.addHistogram(2560.0, 5120.0, bookHistogram1D(28,1,1));
      _difT.addHistogram(5120.0, 10240.0, bookHistogram1D(29,1,1));
      _difT.addHistogram(10240.0, 20480.0, bookHistogram1D(30,1,1));

      ///Broadening
      _difB.addHistogram(80.0, 160.0, bookHistogram1D(31,1,1));
      _difB.addHistogram(160.0, 320.0, bookHistogram1D(32,1,1));
      _difB.addHistogram(320.0, 640.0, bookHistogram1D(33,1,1));
      _difB.addHistogram(640.0, 1280.0, bookHistogram1D(34,1,1));
      _difB.addHistogram(1280.0, 2560.0, bookHistogram1D(35,1,1));
      _difB.addHistogram(2560.0, 5120.0, bookHistogram1D(36,1,1));
      _difB.addHistogram(5120.0, 10240.0, bookHistogram1D(37,1,1));
      _difB.addHistogram(10240.0, 20480.0, bookHistogram1D(38,1,1));
      
      ///Thrust Gamma
      _difTG.addHistogram(80.0, 160.0, bookHistogram1D(39,1,1));
      _difTG.addHistogram(160.0, 320.0, bookHistogram1D(40,1,1));
      _difTG.addHistogram(320.0, 640.0, bookHistogram1D(41,1,1));
      _difTG.addHistogram(640.0, 1280.0, bookHistogram1D(42,1,1));
      _difTG.addHistogram(1280.0, 2560.0, bookHistogram1D(43,1,1));
      _difTG.addHistogram(2560.0, 5120.0, bookHistogram1D(44,1,1));
      _difTG.addHistogram(5120.0, 10240.0, bookHistogram1D(45,1,1));
      _difTG.addHistogram(10240.0, 20480.0, bookHistogram1D(46,1,1));

      ///Broadening Gamma
      _difBG.addHistogram(80.0, 160.0, bookHistogram1D(47,1,1));
      _difBG.addHistogram(160.0, 320.0, bookHistogram1D(48,1,1));
      _difBG.addHistogram(320.0, 640.0, bookHistogram1D(49,1,1));
      _difBG.addHistogram(640.0, 1280.0, bookHistogram1D(50,1,1));
      _difBG.addHistogram(1280.0, 2560.0, bookHistogram1D(51,1,1));
      _difBG.addHistogram(2560.0, 5120.0, bookHistogram1D(52,1,1));
      _difBG.addHistogram(5120.0, 10240.0, bookHistogram1D(53,1,1));
      _difBG.addHistogram(10240.0, 20480.0, bookHistogram1D(54,1,1));

      ///Event Shapes:
      _hist1MinusT = bookHistogram1D("ThrustT",10, 0, 1);
      _hist1MinusB = bookHistogram1D("BroadeningT", 10, 0, 1);
      _histM = bookHistogram1D("InvariantJetMass",10, 0, 1);
      _histC = bookHistogram1D("CParam",10, 0, 1);
      _hist1MinusTG = bookHistogram1D("ThrustG",10, 0, 1);
      _hist1MinusBG = bookHistogram1D("BroadeningG",10, 0, 1);

      ///and profiles:
      ///Mean Event Shape Variables
      _prof1MinusT = bookProfile1D(1,1,1);
      _prof1MinusT2 = bookProfile1D(1,2,1);
      
      _prof1MinusB = bookProfile1D(2,1,1);  
      _prof1MinusB2 = bookProfile1D(2,2,1); 

      _profM = bookProfile1D(3,1,1);
      _profM2 = bookProfile1D(3,2,1);

      _profC = bookProfile1D(4,1,1);
      _profC2 = bookProfile1D(4,2,1);

      _prof1MinusTG = bookProfile1D(5,1,1);
      _prof1MinusTG2 = bookProfile1D(5,2,1);
      
      _prof1MinusBG = bookProfile1D(6,1,1);  
      _prof1MinusBG2 = bookProfile1D(6,2,1); 

      _profQ2vsX = bookProfile1D("Q2vsX",100,0,1);
      _profQ2vsY = bookProfile1D("Q2vsY",100,0,1);
      _profYvsX  = bookProfile1D("YvsX",100,0,1);

    }



    void analyze(const Event& event) {

      double weight = event.weight();

      ///Apply the Projections needed more than once:
      const FinalState& fs = applyProjection<FinalState>(event, "FS");

      const DISKinematics& dk = applyProjection<DISKinematics>(event, "DK");

      const FinalStateBreit& fsb = applyProjection<FinalStateBreit>(event, "FSB");


      ///Get final state information:

      ///Get DIS kinematics:
      ///Q^2: 
      double q2 = dk.Q2();
      ///bjorken x:
      double x = dk.x();
      ///inlasticity:
      double y = dk.y();
      ///W^2
      double w2 = dk.W2();

      ///Apply other Projections:
      ///Invariant Jet Mass
      const InvariantJetMass& invariantjetmassb = 
	applyProjection<InvariantJetMass>(event, "IJM");

      ///C-Parameter
      const Cparam& cparamb = applyProjection<Cparam>(event, "CPM");

      ///Thrust in current hemisphere of Breit frame
      const Thrust& thrustb = applyProjection<Thrust>(event, "Thrust");

      ///Broadening in current hemisphere of Breit Frame
      const Broadening& broadeningb =
	applyProjection<Broadening>(event, "Broadening");

      ///Thrust Gamma in the current hemisphere of the Breit Frame
      const ThrustGamma& thrustgammab =
	applyProjection<ThrustGamma>(event, "ThrustGamma");

      ///Broadening Gamma in the current hemisphere of the Breit Frame
      const BroadeningGamma& broadeninggammab = 
	applyProjection<BroadeningGamma>(event, "BroadeningGamma");
      

      ///before filling, double check we are within Q, x and y limits
      ///- Should have been done in MCGenerator
      if(((q2>80)&&(q2<20480))&&((x>.0024)&&(x<.6))&&((y>.04)&&(y<.9))){

	///DIS Variables
	_histQ -> fill(q2, weight);
	_histX -> fill(x, weight);
	_histY -> fill(y, weight);
	
	///Event Shape Variables
	_difM.fill(q2, invariantjetmassb.invariantjetmass(), weight);

	_difC.fill(q2, cparamb.cparam(), weight);

	_difT.fill(q2, thrustb.thrust(), weight);

	_difB.fill(q2, broadeningb.broadening(), weight);

	_difTG.fill(q2, thrustgammab.thrustgamma(), weight);

	_difBG.fill(q2, broadeninggammab.broadeninggamma(), weight);

	///Event Shapes
	_hist1MinusT -> fill(1.0 - thrustb.thrust(), weight);
	_hist1MinusB -> fill(1.0 - broadeningb.broadening(), weight);
	_histM -> fill(1.0 - invariantjetmassb.invariantjetmass(), weight);
	_histC -> fill(1.0 - cparamb.cparam(), weight);
	_hist1MinusTG -> fill(1.0 - thrustgammab.thrustgamma(), weight);
	_hist1MinusBG -> fill(1.0 -broadeninggammab.broadeninggamma(), weight);

	///Mean Event Shape Variables
	_prof1MinusT -> fill(q2, 1.0 - thrustb.thrust(), weight);
	_prof1MinusT2 -> fill(x, 1.0 - thrustb.thrust(), weight);
	  
	_prof1MinusB -> fill(q2, 1.0 - broadeningb.broadening(), weight);
	_prof1MinusB2 -> fill(x, 1.0 - broadeningb.broadening(), weight);
	
	_profM -> fill(q2, invariantjetmassb.invariantjetmass(), weight);
	_profM2 -> fill(q2, invariantjetmassb.invariantjetmass(), weight);

	_profC -> fill(q2, cparamb.cparam(), weight);
	_profC2 -> fill(x, cparamb.cparam(), weight);

	_prof1MinusTG -> fill(q2, 1.0 - thrustgammab.thrustgamma(), weight);
	_prof1MinusTG2 -> fill(x, 1.0 - thrustgammab.thrustgamma(), weight);
	  
	_prof1MinusBG -> fill(q2, 1.0 - broadeninggammab.broadeninggamma(), weight);
	_prof1MinusBG2 -> fill(x, 1.0 - broadeninggammab.broadeninggamma(), weight);

	///other DIS:

	_profQ2vsX -> fill(x, q2, weight);
	_profQ2vsY -> fill(y, q2, weight);
	_profYvsX  -> fill(x, y, weight);


      }
    }

    
	

    void finalize() {

      foreach (AIDA::IHistogram1D* hist, _difM.getHistograms()) {
	normalize(hist);
      }

      foreach (AIDA::IHistogram1D* hist, _difC.getHistograms()) {
	normalize(hist);
      }

      foreach (AIDA::IHistogram1D* hist, _difT.getHistograms()) {
	normalize(hist);
      }

      foreach (AIDA::IHistogram1D* hist, _difB.getHistograms()) {
	normalize(hist);
      }

      foreach (AIDA::IHistogram1D* hist, _difTG.getHistograms()) {
	normalize(hist);
      }

      foreach (AIDA::IHistogram1D* hist, _difBG.getHistograms()) {
	normalize(hist);
      }

      normalize(_hist1MinusT);
      normalize(_hist1MinusB);
      normalize(_histM);
      normalize(_histC);
      normalize(_hist1MinusTG);
      normalize(_hist1MinusBG);
    }



  private:

    AIDA::IHistogram1D *_histQ;
    AIDA::IHistogram1D *_histX;
    AIDA::IHistogram1D *_histY;

    AIDA::IHistogram1D *_hist1MinusT;
    AIDA::IHistogram1D *_hist1MinusB;
    AIDA::IHistogram1D *_histM;
    AIDA::IHistogram1D *_histC;
    AIDA::IHistogram1D *_hist1MinusTG;
    AIDA::IHistogram1D *_hist1MinusBG;

    AIDA::IProfile1D *_prof1MinusT;
    AIDA::IProfile1D *_prof1MinusT2;
    AIDA::IProfile1D *_prof1MinusB;
    AIDA::IProfile1D *_prof1MinusB2;

    AIDA::IProfile1D *_profM;
    AIDA::IProfile1D *_profM2;
    AIDA::IProfile1D *_profC;
    AIDA::IProfile1D *_profC2;

    AIDA::IProfile1D *_prof1MinusTG;
    AIDA::IProfile1D *_prof1MinusTG2;
    AIDA::IProfile1D *_prof1MinusBG;
    AIDA::IProfile1D *_prof1MinusBG2;

    AIDA::IProfile1D *_profQ2vsX;
    AIDA::IProfile1D *_profQ2vsY;
    AIDA::IProfile1D *_profYvsX;

    

    BinnedHistogram<double> _difM;

    BinnedHistogram<double> _difC;

    BinnedHistogram<double> _difT;

    BinnedHistogram<double> _difB;

    BinnedHistogram<double> _difTG;
   
    BinnedHistogram<double> _difBG;

    
  };



  // This global object acts as a hook for the plugin system
  AnalysisBuilder<ZEUS_2007_S6591329> plugin_ZEUS_2007_S6591329;



}