#include "Rivet/Analysis.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/ChargedFinalState.hh"
#include "Rivet/AnalysisLoader.hh"
#include "Rivet/RivetAIDA.hh"
#include "Rivet/Math/Constants.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"
#include "Rivet/Projections/FinalState.hh"
#include "Rivet/Projections/DISKinematics.hh"
#include "Rivet/Projections/Thrust.hh"
#include "FinalStateBreit.hh"
#include "Broadening.hh"
#include "FinalStateBreitCurrent.hh"
#include "ThrustGamma.hh"
#include "InvariantJetMass.hh"
#include "Cparam.hh"

namespace Rivet {

  ///@brief ZEUS Event shapes in deep inelastic scattering at HERA
  ///@author Alexander Radovic
  ///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);
    }

   
    void init() {
      /// Add the projections used from Rivet library
      const FinalState fs;
      addProjection(fs,"FS");

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

      const Thrust thrust(fs);
      addProjection(thrust, "Thrust");
    
      ///Add the additional projections created for this analysis
      const FinalStateBreitCurrent fsb= FinalStateBreitCurrent(dk,1,.25);
      addProjection(fsb, "FSB");

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

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

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

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

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

      ///Book histograms for the mean variables
      
      _hist1MinusT2 = bookHistogram1D("Thrust2",10, 0, 1);
      _hist1MinusB = bookHistogram1D("Broadening",10, 0, 1);
      _hist1M = bookHistogram1D("InvariantJetMass",10, 0, 1);
      _hist1C = bookHistogram1D("CParam",10, 0, 1);


      ///Book the histograms for differential distributions 
      ///using the referance data:

      //thust gamma histos
      _difTr1G= bookHistogram1D(39,1,1);
      _difTr2G= bookHistogram1D(40,1,1);
      _difTr3G= bookHistogram1D(41,1,1);
      _difTr4G= bookHistogram1D(42,1,1);
      _difTr5G= bookHistogram1D(43,1,1);
      _difTr6G= bookHistogram1D(44,1,1);
      _difTr7G= bookHistogram1D(45,1,1);
      _difTr8G= bookHistogram1D(46,1,1);

      //broadening gamma histos
      _difBr1G= bookHistogram1D(47,1,1);
      _difBr2G= bookHistogram1D(48,1,1);
      _difBr3G= bookHistogram1D(49,1,1);
      _difBr4G= bookHistogram1D(50,1,1);
      _difBr5G= bookHistogram1D(51,1,1);
      _difBr6G= bookHistogram1D(52,1,1);
      _difBr7G= bookHistogram1D(53,1,1);
      _difBr8G= bookHistogram1D(54,1,1);
   
      //thrust histos
      _difTr1= bookHistogram1D(23,1,1);
      _difTr2= bookHistogram1D(24,1,1);
      _difTr3= bookHistogram1D(25,1,1);
      _difTr4= bookHistogram1D(26,1,1);
      _difTr5= bookHistogram1D(27,1,1);
      _difTr6= bookHistogram1D(28,1,1);
      _difTr7= bookHistogram1D(29,1,1);
      _difTr8= bookHistogram1D(30,1,1);

      //Invariant Jet Mass histos
      _difM1= bookHistogram1D(7,1,1);
      _difM2= bookHistogram1D(8,1,1);
      _difM3= bookHistogram1D(9,1,1);
      _difM4= bookHistogram1D(10,1,1);
      _difM5= bookHistogram1D(11,1,1);
      _difM6= bookHistogram1D(12,1,1);
      _difM7= bookHistogram1D(13,1,1);
      _difM8= bookHistogram1D(14,1,1);

      //C Parameter histos
      _difC1= bookHistogram1D(15,1,1);
      _difC2= bookHistogram1D(16,1,1);
      _difC3= bookHistogram1D(17,1,1);
      _difC4= bookHistogram1D(18,1,1);
      _difC5= bookHistogram1D(19,1,1);
      _difC6= bookHistogram1D(20,1,1);
      _difC7= bookHistogram1D(21,1,1);
      _difC8= bookHistogram1D(22,1,1);

      //histograms of the broadening
      _difBr1= bookHistogram1D(31,1,1);
      _difBr2= bookHistogram1D(32,1,1);
      _difBr3= bookHistogram1D(33,1,1);
      _difBr4= bookHistogram1D(34,1,1);
      _difBr5= bookHistogram1D(35,1,1);
      _difBr6= bookHistogram1D(36,1,1);
      _difBr7= bookHistogram1D(37,1,1);
      _difBr8= bookHistogram1D(38,1,1);
      

      ///and for the DIS variables

      _histQ = bookHistogram1D("Q2",100, 0, 20480);
      _histX = bookHistogram1D("X",100, 0,1);
      _histY = bookHistogram1D("Y",100, 0,1);

      ///book Profile1D for variables that arn't specific to one q2 bin
      ///allows us to plot against other variables

      _profQ2vsX  = bookProfile1D("Q2vsX",100,0,1);
      _profQ2vsY  = bookProfile1D("Q2vsY",100,0,1);
      _profYvsX   = bookProfile1D("YvsX",100,0,1);
      
      _profQ2vsMT = bookProfile1D(1,1,1);
      _profQ2vsB  = bookProfile1D(2,1,1);
      _profQ2vsM2 = bookProfile1D(3,1,1);
      _profQ2vsC  = bookProfile1D(4,1,1);
      
    }
    
    void analyze(const Event& event) {

      /// Apply projections which will be re-used
      const FinalState& fs = applyProjection<FinalState>(event, "FS");
      const DISKinematics& dk = applyProjection<DISKinematics>(event, "DK");
      
      const FinalStateBreitCurrent& fsb = applyProjection<FinalStateBreitCurrent>(event, "FSB");
     
      ///get the weight for each event:     
      const double weight = event.weight();
      
      //Get final state information:
      
      //Get DIS kinematics:
      //Q^2: 
      double q2 = dk.Q2();
      //bjorken:
      double x = dk.x();
      //inlasticity:
      double y = dk.y();
      //W^2
      double w2 = dk.W2();
      
          
           
      ///Thrust in the current hemisphere of the breit frame:
      const Thrust& thrustb = applyProjection<Thrust>(event, "BThrust");
      _hist1MinusT2->fill(1.0 - thrustb.thrust(), weight);

      ///Thrust along the Gamma axis in breit frame
      const ThrustGamma& thrustgamma = applyProjection<ThrustGamma>(event, "ThrustGamma");

      ///Broadening in current hemisphere of breit frame
      const Broadening& broadeningb = applyProjection<Broadening>(event, "Broadening");
      _hist1MinusB->fill(broadeningb.broadening(), weight);

      ///Apply Invariant Jet Mass Projection and fill mean value histogram
      const InvariantJetMass& invariantjetmassb = applyProjection<InvariantJetMass>(event, "IJM");
      _hist1M->fill(invariantjetmassb.invariantjetmass(), weight);

      ///Apply C-paramater Projection and fill histogram
      const Cparam& cparamb = applyProjection<Cparam>(event, "CPM");
      _hist1C->fill(cparamb.cparam(), weight);


      ///check within q2, x & y limits:
      ///These cuts should be made in the MC Generator to avoid waste data

      if((q2 > 80) && (q2 < 20480) && ((x > .0024) && (x < .6)) && ((y > .04) && (y < .9))){
	  _histQ->fill(q2, weight);
	  _histX->fill(x, weight);
	  _histY->fill(y, weight);

	  /// work out average thrusts etc for differant q^2 bins:
	  if (q2 < 160){
	    _difTr1->fill(thrustb.thrust(), event.weight());
	    _difBr1->fill(broadeningb.broadening(), event.weight());
	    _difTr1G->fill(thrustb.thrust(), event.weight());
	    _difBr1G->fill(broadeningb.broadening(), event.weight());
	    _difM1->fill(invariantjetmassb.invariantjetmass(), event.weight());
	    _difC1->fill(cparamb.cparam(), event.weight());
	  }
	  else if (q2 < 320){
	    _difTr2->fill(thrustb.thrust(), event.weight());
	    _difBr2->fill(broadeningb.broadening(), event.weight());	
	    _difTr2G->fill(thrustb.thrust(), event.weight());
	    _difBr2G->fill(broadeningb.broadening(), event.weight());
	    _difM2->fill(invariantjetmassb.invariantjetmass(), event.weight());
	    _difC2->fill(cparamb.cparam(), event.weight());
	  }
	  else if (q2 < 640){	
	    _difTr3->fill(thrustb.thrust(), event.weight());
	    _difBr3->fill(broadeningb.broadening(), event.weight());
	    _difTr3G->fill(thrustb.thrust(), event.weight());
	    _difBr3G->fill(broadeningb.broadening(), event.weight());
	    _difM3->fill(invariantjetmassb.invariantjetmass(), event.weight());
	    _difC3->fill(cparamb.cparam(), event.weight());
	  }
	  else if (q2 < 1280){
	    _difTr4->fill(thrustb.thrust(), event.weight());
	    _difBr4->fill(broadeningb.broadening(), event.weight());
	    _difTr4G->fill(thrustb.thrust(), event.weight());
	    _difBr4G->fill(broadeningb.broadening(), event.weight());
	    _difM4->fill(invariantjetmassb.invariantjetmass(), event.weight());
	    _difC4->fill(cparamb.cparam(), event.weight());
	  }
	  else if (q2 < 2560){
	    _difTr5->fill(thrustb.thrust(), event.weight());
	    _difBr5->fill(broadeningb.broadening(), event.weight());
	    _difTr5G->fill(thrustb.thrust(), event.weight());
	    _difBr5G->fill(broadeningb.broadening(), event.weight());
	    _difM5->fill(invariantjetmassb.invariantjetmass(), event.weight());
	    _difC5->fill(cparamb.cparam(), event.weight());
	  }
	  else if (q2 < 5120){
	    _difTr6->fill(thrustb.thrust(), event.weight());
	    _difBr6->fill(broadeningb.broadening(), event.weight());
	    _difTr6G->fill(thrustb.thrust(), event.weight());
	    _difBr6G->fill(broadeningb.broadening(), event.weight());
	    _difM6->fill(invariantjetmassb.invariantjetmass(), event.weight());
	    _difC6->fill(cparamb.cparam(), event.weight());
	  }
	  else if (q2 < 10240){
	    _difTr7->fill(thrustb.thrust(), event.weight());
	    _difBr7->fill(broadeningb.broadening(), event.weight());
	    _difTr7G->fill(thrustb.thrust(), event.weight());
	    _difBr7G->fill(broadeningb.broadening(), event.weight());
	    _difM7->fill(invariantjetmassb.invariantjetmass(), event.weight());
	    _difC7->fill(cparamb.cparam(), event.weight());
	  }
	  else if (q2 < 20480){
	    _difTr8->fill(thrustb.thrust(), event.weight());
	    _difBr8->fill(broadeningb.broadening(), event.weight());
	    _difTr8G->fill(thrustb.thrust(), event.weight());
	    _difBr8G->fill(broadeningb.broadening(), event.weight());
	    _difM8->fill(invariantjetmassb.invariantjetmass(), event.weight());
	    _difC8->fill(cparamb.cparam(), event.weight());
	  }

	  ///put data to all the plots that arent for a specific Q2 bin:
	  _profQ2vsMT->fill(q2,1.0-thrustb.thrust(), event.weight() );
	  _profQ2vsB->fill(q2,broadeningb.broadening(), event.weight() );
	  _profQ2vsC->fill(q2,cparamb.cparam(), event.weight() );
	  _profQ2vsM2->fill(q2,invariantjetmassb.invariantjetmass(), event.weight() );

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

    void finalize() {

      ///normalize all my histograms:
       normalize(_hist1MinusT);
       normalize(_hist1MinusT2);
       normalize(_hist1MinusB);
       normalize(_hist1M);
       normalize(_hist1C);

      
       normalize(_difM1);
       normalize(_difM2);
       normalize(_difM3);
       normalize(_difM4);
       normalize(_difM5);
       normalize(_difM6);
       normalize(_difM7);
       normalize(_difM8);

       normalize(_difC1);
       normalize(_difC2);
       normalize(_difC3);
       normalize(_difC4);
       normalize(_difC5);
       normalize(_difC6);
       normalize(_difC7);
       normalize(_difC8);

       normalize(_difTr1);
       normalize(_difTr2);
       normalize(_difTr3);
       normalize(_difTr4);
       normalize(_difTr5);
       normalize(_difTr6);
       normalize(_difTr7);
       normalize(_difTr8);

       normalize(_difTr1G);
       normalize(_difTr2G);
       normalize(_difTr3G);
       normalize(_difTr4G);
       normalize(_difTr5G);
       normalize(_difTr6G);
       normalize(_difTr7G);
       normalize(_difTr8G);

       normalize(_difBr1);
       normalize(_difBr2);
       normalize(_difBr3);
       normalize(_difBr4);
       normalize(_difBr5);
       normalize(_difBr6);
       normalize(_difBr7);
       normalize(_difBr8);

       normalize(_difBr1G);
       normalize(_difBr2G);
       normalize(_difBr3G);
       normalize(_difBr4G);
       normalize(_difBr5G);
       normalize(_difBr6G);
       normalize(_difBr7G);
       normalize(_difBr8G);
  
       
  }


    private: 
    
    AIDA::IHistogram1D *_hist1MinusT;
    AIDA::IHistogram1D *_hist1MinusT2;
    AIDA::IHistogram1D *_hist1MinusB;
    AIDA::IHistogram1D *_hist1M;
    AIDA::IHistogram1D *_hist1C;
    AIDA::IHistogram1D *_histQ;
    AIDA::IHistogram1D *_histX;
    AIDA::IHistogram1D *_histY;


    AIDA::IHistogram1D *_difTr1;
    AIDA::IHistogram1D *_difTr2;
    AIDA::IHistogram1D *_difTr3;
    AIDA::IHistogram1D *_difTr4;
    AIDA::IHistogram1D *_difTr5;
    AIDA::IHistogram1D *_difTr6;
    AIDA::IHistogram1D *_difTr7;
    AIDA::IHistogram1D *_difTr8;


    AIDA::IHistogram1D *_difM1;
    AIDA::IHistogram1D *_difM2;
    AIDA::IHistogram1D *_difM3;
    AIDA::IHistogram1D *_difM4;
    AIDA::IHistogram1D *_difM5;
    AIDA::IHistogram1D *_difM6;
    AIDA::IHistogram1D *_difM7;
    AIDA::IHistogram1D *_difM8;

    AIDA::IHistogram1D *_difC1;
    AIDA::IHistogram1D *_difC2;
    AIDA::IHistogram1D *_difC3;
    AIDA::IHistogram1D *_difC4;
    AIDA::IHistogram1D *_difC5;
    AIDA::IHistogram1D *_difC6;
    AIDA::IHistogram1D *_difC7;
    AIDA::IHistogram1D *_difC8;

    AIDA::IHistogram1D *_difTr1G;
    AIDA::IHistogram1D *_difTr2G;
    AIDA::IHistogram1D *_difTr3G;
    AIDA::IHistogram1D *_difTr4G;
    AIDA::IHistogram1D *_difTr5G;
    AIDA::IHistogram1D *_difTr6G;
    AIDA::IHistogram1D *_difTr7G;
    AIDA::IHistogram1D *_difTr8G;

    AIDA::IHistogram1D *_difBr1;
    AIDA::IHistogram1D *_difBr2;
    AIDA::IHistogram1D *_difBr3;
    AIDA::IHistogram1D *_difBr4;
    AIDA::IHistogram1D *_difBr5;
    AIDA::IHistogram1D *_difBr6;
    AIDA::IHistogram1D *_difBr7;
    AIDA::IHistogram1D *_difBr8;

    AIDA::IHistogram1D *_difBr1G;
    AIDA::IHistogram1D *_difBr2G;
    AIDA::IHistogram1D *_difBr3G;
    AIDA::IHistogram1D *_difBr4G;
    AIDA::IHistogram1D *_difBr5G;
    AIDA::IHistogram1D *_difBr6G;
    AIDA::IHistogram1D *_difBr7G;
    AIDA::IHistogram1D *_difBr8G;

    AIDA::IProfile1D *_profQ2vsMT;
    AIDA::IProfile1D *_profQ2vsB;
    AIDA::IProfile1D *_profQ2vsM2;
    AIDA::IProfile1D *_profQ2vsC;
    AIDA::IProfile1D *_profQ2vsX;
    AIDA::IProfile1D *_profQ2vsY;
    AIDA::IProfile1D *_profYvsX;


  };

  AnalysisBuilder<ZEUS_2007_S6591329> plugin_ZEUS_2007_S6591329;

}