// -*- C++ -*-
//@author Peter Wijeratne <paw@hep.ucl.ac.uk>
//@author Robindra Prabhu <prabhu@cern.ch>
#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/VetoedFinalState.hh"
#include "Rivet/Projections/FastJets.hh"
#include "Rivet/Tools/ParticleIdUtils.hh"

#include "math.h"


namespace Rivet{
  
  ///A very basic analysis sensitive to ET flow in minbias and dijet events
  class ATLAS_2012_I1183818 : public Analysis {

  public:
    
    ATLAS_2012_I1183818()
      : Analysis("ATLAS_2012_I1183818")
    {}
    
    ///////////////////////////book histograms and initialise projections

  public:
    
    void init(){

      const FinalState cnfs(-4.8,4.8, 0*MeV);
      const ChargedFinalState cfs(-2.5, 2.5, 250*MeV);
      addProjection(cnfs, "FS");
      addProjection(cfs, "CFS");

      //handle on jets
      const FastJets jetsAntiKt4(cnfs, FastJets::ANTIKT, 0.4);
      addProjection(jetsAntiKt4, "AntiKt4Jets");

      // ------- MINBIAS HISTOGRAMS --------
      //
      ///MB event counter                                                        
      m_chargedEvents = 0.0;

      _histETflowEtaNorm = bookHistogram1D(1,1,1);
      _histSumETbin1 = bookHistogram1D(3,1,1);
      _histSumETbin2 = bookHistogram1D(4,1,1);
      _histSumETbin3 = bookHistogram1D(5,1,1);
      _histSumETbin4 = bookHistogram1D(6,1,1);
      _histSumETbin5 = bookHistogram1D(7,1,1);
      _histSumETbin6 = bookHistogram1D(8,1,1);
      


      // ------- DIJET HISTOGRAMS --------
      //
      // Do dijet analysis?
      doDijetAnalysis = true;
      // Dijet event counter                                                                                                     
      m_events_dijets = 0.0;

      // Optional helper histograms to hold the number of events that 1) were processed, 2) passed the MB selection, 3) passed the dijet selection
      // ---> include if required for normalising multiple output jobs

      //_histEventsProcessed = bookHistogram1D("Counter_EventsProcessed",2, 0.5, 2.5 );
      //_histEventsPassedMB  = bookHistogram1D("Counter_EventsPassedMB",2, 0.5, 2.5 );
      //_histEventsPassedUE  = bookHistogram1D("Counter_EventsPassedUE",2, 0.5, 2.5 );


      //sumET
      _histETflowTransverseAbsEtaNorm = bookHistogram1D(2,1,1);

      _histSumETTransverseAbsEtaNorm_VarBins_0p0_0p8 = bookHistogram1D(9,1,1);
      _histSumETTransverseAbsEtaNorm_VarBins_0p8_1p6 = bookHistogram1D(10,1,1);
      _histSumETTransverseAbsEtaNorm_VarBins_1p6_2p4 = bookHistogram1D(11,1,1);
      _histSumETTransverseAbsEtaNorm_VarBins_2p4_3p2 = bookHistogram1D(12,1,1);
      _histSumETTransverseAbsEtaNorm_VarBins_3p2_4p0 = bookHistogram1D(13,1,1);
      _histSumETTransverseAbsEtaNorm_VarBins_4p0_4p8 = bookHistogram1D(14,1,1);


    }///end of init

    ///////////////////////////loop over events 

    void analyze(const Event& event){
      
      const double weight = event.weight();

      //_histEventsProcessed->fill(1, weight); // uncomment this to include histogram holding number of events processed

      const FinalState& cfs = applyProjection<FinalState>(event, "CFS");

      bool isCharged = false;
      //event selection: > 2 charged particles with pT > 250.MeV and |eta| < 2.5
      if(cfs.size() > 1){
	isCharged = true;
	m_chargedEvents++;
	//_histEventsPassedMB->fill(1,  weight); // uncomment this to include histogram holding number of events passing the MB selection
      }
      
      const FinalState& cnfs = applyProjection<FinalState>(event, "FS");
      //      getLog() << Log::DEBUG << "Total multiplicity = " << cnfs.size() << endl;
      
      //get jets
      const FastJets& jetsAntiKt4 = applyProjection<FastJets>(event, "AntiKt4Jets");
      const Jets& jAkt4 = jetsAntiKt4.jetsByPt(20.0*GeV);
      
      //numAntiKt4Jets += jAkt4.size();
      
      //initialise sumET variables
      double sumETbin1 = 0;
      double sumETbin2 = 0;
      double sumETbin3 = 0;
      double sumETbin4 = 0;
      double sumETbin5 = 0;
      double sumETbin6 = 0;
      //initialise dijet sumET variables
      m_trans_sumET_0p0_0p8 = 0.;
      m_trans_sumET_0p8_1p6 = 0.;
      m_trans_sumET_1p6_2p4 = 0.;
      m_trans_sumET_2p4_3p2 = 0.;
      m_trans_sumET_3p2_4p0 = 0.;
      m_trans_sumET_4p0_4p8 = 0.;
      
      
      // if(passes event selection)
      if(isCharged){
	
	foreach( const Particle& p, cnfs.particles() ){

	  double pp = p.momentum().p().mod();
	  
	  ///enforce truth selection representing detected particle sensitivity

	  if(pp < 0.5*GeV && PID::threeCharge(p.pdgId()) != 0) continue;
  
	  if(pp < 0.2*GeV && PID::threeCharge(p.pdgId()) == 0) continue;

	  ///

	  ///calculate variables
	  double ET = (p.momentum().E()/GeV)*sin( 2.*atan(exp(-1.*p.momentum().eta())));
	  double feta = std::fabs(p.momentum().eta());

	  ///fill histograms
	  _histETflowEtaNorm->fill(feta, weight*ET); 
	  
	  if(feta <= 0.8) sumETbin1 += ET; 
	  if(feta > 0.8 && feta <= 1.6) sumETbin2 += ET;
	  if(feta > 1.6 && feta <= 2.4) sumETbin3 += ET; 
	  if(feta > 2.4 && feta <= 3.2) sumETbin4 += ET; 
	  if(feta > 3.2 && feta <= 4.0) sumETbin5 += ET; 
	  if(feta > 4.0) sumETbin6 += ET; 
	  
	}//end of foreach
	
	_histSumETbin1->fill(sumETbin1,weight);
	_histSumETbin2->fill(sumETbin2,weight);
	_histSumETbin3->fill(sumETbin3,weight);
	_histSumETbin4->fill(sumETbin4,weight);
	_histSumETbin5->fill(sumETbin5,weight);
	_histSumETbin6->fill(sumETbin6,weight);
      }

      if( ! doDijetAnalysis ) return;

      // --- do dijet analysis ---

      //require at least two jets
      if( jAkt4.size() < 2 ) { return; } 

      //require two leading jets to pass ET cuts
      if( jAkt4[0].momentum().Et() < 20.*GeV) { return; }
      if( jAkt4[1].momentum().Et() < 20.*GeV) { return; }

      //require leading jets to be central
      if( fabs( jAkt4[0].momentum().eta() ) > 2.5 ) { return; }
      if( fabs( jAkt4[1].momentum().eta() ) > 2.5 ) { return; }

      //require back-to-back topology
      float dphi = deltaPhi( jAkt4[0].momentum(), jAkt4[1].momentum() );
      if( fabs(dphi) < 2.5 ) { return; }

      //require ET-balance
      float etbalance = ( jAkt4[1].momentum().Et() / jAkt4[0].momentum().Et() );
      if( etbalance < 0.5 ) { return; }

      //found an event that satisfies dijet selection, now fill histograms...
      fillUEHistograms( jAkt4, cnfs, weight);
      fillUESumETHistograms(weight);

    }///end of analyze
    
    ///////////////////////////finalize histograms - scale, divide etc.
    
    void finalize(){

      /// automatically scales by binwidth and event weight = 1
      /// several scale factors here:
      /// 1. nEvents (m_chargedEvents)
      /// 2. phase-space (2*M_PI)
      /// 3. double binning due to symmetrisation (2)
      scale( _histETflowEtaNorm, ((double)1/(m_chargedEvents*4.*M_PI)) );      
      scale( _histSumETbin1,     ((double)1/(m_chargedEvents)) );
      scale( _histSumETbin2,     ((double)1/(m_chargedEvents)) );
      scale( _histSumETbin3,     ((double)1/(m_chargedEvents)) );
      scale( _histSumETbin4,     ((double)1/(m_chargedEvents)) );
      scale( _histSumETbin5,     ((double)1/(m_chargedEvents)) );
      scale( _histSumETbin6,     ((double)1/(m_chargedEvents)) );

      //Dijet analysis
      if( ! doDijetAnalysis ) return;

      // Dijet scale factors:
      //1. number of events passing dijet selection
      //2. phase-space: 1. / 2/3*M_PI
      //3. double binning due to symmetrisation in |eta| plots : 1/2 

      getLog() << Log::INFO <<"Number of events passing dijet selection:  "<< m_events_dijets <<endl;

      double scaleFactor = (1. / m_events_dijets) * (1. / ( (2./3.)*M_PI ) );
    
      //Transverse region
      scale( _histETflowTransverseAbsEtaNorm, scaleFactor*0.5 );

      scale( _histSumETTransverseAbsEtaNorm_VarBins_0p0_0p8, (1. / m_events_dijets) );
      scale( _histSumETTransverseAbsEtaNorm_VarBins_0p8_1p6, (1. / m_events_dijets) );
      scale( _histSumETTransverseAbsEtaNorm_VarBins_1p6_2p4, (1. / m_events_dijets) );
      scale( _histSumETTransverseAbsEtaNorm_VarBins_2p4_3p2, (1. / m_events_dijets) );
      scale( _histSumETTransverseAbsEtaNorm_VarBins_3p2_4p0, (1. / m_events_dijets) );
      scale( _histSumETTransverseAbsEtaNorm_VarBins_4p0_4p8, (1. / m_events_dijets) );


    }

    void fillUEHistograms(const Jets& hardScale, const FinalState& softScale, double weight) {


      m_events_dijets++;
      //_histEventsPassedUE ->fill( 1, weight); // uncomment this to include histogram holding number of events passing the dijet selection

      //loop over all particles and check their relation to hardscale object
      foreach( const Particle& particle, softScale.particles() ){

	///enforce truth selection representing detected particle sensitivity                              
	double pp = particle.momentum().p().mod();
	// require charged particles to have |p|>500 MeV
	if(pp < 0.5*GeV && PID::threeCharge(particle.pdgId()) != 0) continue;
	// require neutral particles to have |p|>200 MeV
	if(pp < 0.2*GeV && PID::threeCharge(particle.pdgId()) == 0) continue;
	///                                                                         

	///calculate variables
	double dPhi = deltaPhi( hardScale[0], particle.momentum() );
	double ET   = (particle.momentum().E())*sin( 2.*atan(exp(-1.*particle.momentum().eta())));
	ET          = ET*GeV;
	double Eta  = particle.momentum().eta();


	//Transverse region
	if( ( fabs(dPhi) > (1./3.)*M_PI ) && ( fabs(dPhi) < (2./3.)*M_PI ) ){
	  
          _histETflowTransverseAbsEtaNorm->fill( fabs(Eta), weight*ET );

	  if(fabs(Eta) >= 0.0 && fabs(Eta) <  0.8){ m_trans_sumET_0p0_0p8 += ET; }
	  if(fabs(Eta) >= 0.8 && fabs(Eta) <  1.6){ m_trans_sumET_0p8_1p6 += ET; }
	  if(fabs(Eta) >= 1.6 && fabs(Eta) <  2.4){ m_trans_sumET_1p6_2p4 += ET; }
	  if(fabs(Eta) >= 2.4 && fabs(Eta) <  3.2){ m_trans_sumET_2p4_3p2 += ET; }
	  if(fabs(Eta) >= 3.2 && fabs(Eta) <  4.0){ m_trans_sumET_3p2_4p0 += ET; }
	  if(fabs(Eta) >= 4.0 && fabs(Eta) <= 4.8){ m_trans_sumET_4p0_4p8 += ET; }

	}
	

      }//end loop over softScale particles

    }

    void fillUESumETHistograms(double weight){

      //Transverse region
      _histSumETTransverseAbsEtaNorm_VarBins_0p0_0p8->fill( m_trans_sumET_0p0_0p8 , weight);
      _histSumETTransverseAbsEtaNorm_VarBins_0p8_1p6->fill( m_trans_sumET_0p8_1p6 , weight);
      _histSumETTransverseAbsEtaNorm_VarBins_1p6_2p4->fill( m_trans_sumET_1p6_2p4 , weight);
      _histSumETTransverseAbsEtaNorm_VarBins_2p4_3p2->fill( m_trans_sumET_2p4_3p2 , weight);
      _histSumETTransverseAbsEtaNorm_VarBins_3p2_4p0->fill( m_trans_sumET_3p2_4p0 , weight);
      _histSumETTransverseAbsEtaNorm_VarBins_4p0_4p8->fill( m_trans_sumET_4p0_4p8 , weight);

    }




  private:

    //Minbias-analysis: variable + histograms
    double m_chargedEvents;

    AIDA::IHistogram1D* _histETflowEtaNorm;
    AIDA::IHistogram1D* _histSumETbin1;
    AIDA::IHistogram1D* _histSumETbin2;
    AIDA::IHistogram1D* _histSumETbin3;
    AIDA::IHistogram1D* _histSumETbin4;
    AIDA::IHistogram1D* _histSumETbin5;
    AIDA::IHistogram1D* _histSumETbin6;


    //Dijet-analysis: variables + histograms
    bool doDijetAnalysis;
    double m_events_dijets;
    AIDA::IHistogram1D *_histEventsProcessed;
    AIDA::IHistogram1D *_histEventsPassedMB; 
    AIDA::IHistogram1D *_histEventsPassedUE;

    //Transverse region
    double m_trans_sumET_0p0_0p8;
    double m_trans_sumET_0p8_1p6;
    double m_trans_sumET_1p6_2p4;
    double m_trans_sumET_2p4_3p2;
    double m_trans_sumET_3p2_4p0;
    double m_trans_sumET_4p0_4p8;
    
    AIDA::IHistogram1D *_histETflowTransverseAbsEtaNorm;
    AIDA::IHistogram1D *_histSumETTransverseAbsEtaNorm_VarBins_0p0_0p8;
    AIDA::IHistogram1D *_histSumETTransverseAbsEtaNorm_VarBins_0p8_1p6;
    AIDA::IHistogram1D *_histSumETTransverseAbsEtaNorm_VarBins_1p6_2p4;
    AIDA::IHistogram1D *_histSumETTransverseAbsEtaNorm_VarBins_2p4_3p2;
    AIDA::IHistogram1D *_histSumETTransverseAbsEtaNorm_VarBins_3p2_4p0;
    AIDA::IHistogram1D *_histSumETTransverseAbsEtaNorm_VarBins_4p0_4p8;
    
  };///end of class ATLAS_2012_I1183818

  DECLARE_RIVET_PLUGIN(ATLAS_2012_I1183818);
  
} ///end of namespace Rivet