#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/AlgFactory.h"
#include "GaudiKernel/IToolSvc.h"
#include "GaudiKernel/GaudiException.h"
#include "GaudiKernel/Property.h"

#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/PropertyMgr.h"
#include <cmath>

#include "GaudiKernel/ITHistSvc.h"
#include "TTree.h"

#include "StoreGate/StoreGateSvc.h"
#include "StoreGate/DataHandle.h"

#include "TH1.h"
#include "TH2.h"

/// the Electron
#include "egammaEvent/Electron.h"
#include "egammaEvent/EMShower.h"
#include "egammaEvent/EMTrackMatch.h"
#include "egammaEvent/egammaPIDdefs.h"
#include "egammaEvent/egammaParamDefs.h"
#include "egammaEvent/egDetailContainer.h"
#include "egammaEvent/egDetail.h"

/// Constituent navigation
#include "Navigation/NavigationToken.h"

/// common implementation of all particles
#include "ParticleEvent/ParticleBaseContainer.h"

/// the composite particle
#include "ParticleEvent/CompositeParticle.h"

/// particle jets
#include "JetEvent/JetCollection.h"

/// analysis tools
#include "AnalysisUtils/AnalysisCombination.h"
#include "UserAnalysisUtils/UserAnalysisSelectionTool.h"

#include "TrigDecisionTool/ChainGroup.h"
#include "TrigDecisionTool/FeatureContainer.h"
#include "TrigDecisionTool/Feature.h"
#include "TrigSteeringEvent/TrigRoiDescriptor.h"
#include "TVector3.h"

#include "VxVertex/VxContainer.h"

/// the header file
#include "AnalysisExamples/MCProduction.h"

#include "McParticleEvent/TruthParticleContainer.h"
#include "TGraph.h"

//#include "AnalysisExamples/checkOQ.h"
//#include "AnalysisExamples/checkOQ.C"
//#include "AnalysisExamples/IsEMPrime.C"

#include <stdint.h>
#include <algorithm>
#include <math.h>
#include <functional>
#include <fstream>
#include <iostream>

using namespace Analysis;
//typedef CompositeParticle::IParticleLink_t IPLink_t;


/////////////////////////////////////////////////////////////////////////////////////
/// Constructor

MCProduction::MCProduction(const std::string& name,
  ISvcLocator* pSvcLocator) : Algorithm(name, pSvcLocator),
  m_analysisTools( "AnalysisTools", this )
{
  /// switches to control the analysis through job options :: these are the default
  /// to changed in the job options

  /// The Electron AOD container name & selection cuts
  declareProperty( "AnalysisTools", m_analysisTools );
  declareProperty("ElectronContainer", m_electronContainerName = "ElectronAODCollection");
  declareProperty("TrackMatchContainer", m_trkMatchContainerName="egDetailContainer");
  //  declareProperty("McParticleContainer", m_truthParticleContainerName = "SpclMC");
  declareProperty("ElectronPtCut", m_ptElecCut = 20.0*GeV);
  declareProperty("ElectronMassCut", m_MassElecCut = 50.0*GeV);
  declareProperty("ElectronEtaCut", m_etaElecCut = 2.47);
  declareProperty("ElectronTriggerFlag", m_ElectronTriggerFlag = "L1_EM10");
  declareProperty("MCorData", m_MCorData = "Data");
  declareProperty("LatestOQ", m_LatestOQ = 165589);
  declareProperty("Scale", m_Scale = 1);
 }

//
/////////////////////////////////////////////////////////////////////////////////////
/// Destructor - check up memory allocation
/// delete any memory allocation on the heap

MCProduction::~MCProduction() {}

////////////////////////////////////////////////////////////////////////////////////
/// Initialize
/// initialize StoreGate
/// get a handle on the analysis tools
/// book histograms

StatusCode MCProduction::initialize() {

  MsgStream mLog( messageService(), name() );

  mLog << MSG::INFO
       << "Initializing MCProduction"
       << endreq;

  StatusCode sc = service("StoreGateSvc", m_storeGate);
  if (sc.isFailure()) {
     mLog << MSG::ERROR
          << "Unable to retrieve pointer to StoreGateSvc"
          << endreq;
     return sc;
  }
  
  /// get a handle on the analysis tools
  sc = m_analysisTools.retrieve();
  if ( sc.isFailure() ) {
      mLog << MSG::ERROR << "Can't get handle on analysis tools" << endreq;
      return sc;
  }

  /// Return a pointer to THistSvc
  sc = service("THistSvc", m_thistSvc);
  if(sc.isFailure() ){
    mLog << MSG::ERROR
         << "Unable to retrieve pointer to THistSvc"
         << endreq;
   return sc;
  }

  // retrieve the TrigDecisionTool
  sc = m_trigDec.retrieve();
  if ( sc.isFailure() ) {
   mLog << MSG::ERROR << "Could not retrieve TrigDecisionTool!" << endreq;
   return StatusCode::FAILURE;
  }

  ///////////////////////////////////////////BOOKING HISTOGRAMS HERE//////////////////////////////////////

  // CutFlow

  h_CutFlow = new TH1F("CutFlow","Cut Flow",17,0,17);
  sc = m_thistSvc->regHist("/AANT/Electron/CutFlow", h_CutFlow);

  /// Electron histogram booking

  // Pair Booking Histograms

  h_Pair_Pt = new TH1F("Pair_Pt","Pair Pt Spectrum",400,0,2000);
  sc = m_thistSvc->regHist("/AANT/Electron/Pair_Pt", h_Pair_Pt);

  h_Pair_Et = new TH1F("Pair_Et","Pair Et Spectrum",400,0,2000);
  sc = m_thistSvc->regHist("/AANT/Electron/Pair_Et", h_Pair_Et);

  h_Pair_Eta = new TH1F("Pair_Eta","Pair Eta Spectrum",60,-3.0,3.0);
  sc = m_thistSvc->regHist("/AANT/Electron/Pair_Eta", h_Pair_Eta);

  h_Pair_InvMass = new TH1F("Pair_InvMass","Pair Invariant Mass Spectrum",800,0,4000);
  sc = m_thistSvc->regHist("/AANT/Electron/Pair_InvMass", h_Pair_InvMass);

  if( sc.isFailure() ){
      mLog << MSG::ERROR << "ROOT Hist aod_zmm_mass_hist registration failed" << endreq;
  }
  return StatusCode::SUCCESS;
}		 

///////////////////////////////////////////////////////////////////////////////////
/// Finalize - delete any memory allocation from the heap

StatusCode MCProduction::finalize() {
  MsgStream mLog( messageService(), name() );
  
  return StatusCode::SUCCESS;

}

//////////////////////////////////////////////////////////////////////////////////
/// Execute - called by the event loop on event by event

StatusCode MCProduction::execute() {

  MsgStream mLog( messageService(), name() );

  mLog << MSG::DEBUG << "execute()" << endreq;

  StatusCode sc = StatusCode::SUCCESS;

  /// do the Z'->ee reconstruction on AOD
  sc = mcproduction();
  if ( sc.isFailure() ) {
    mLog << MSG::ERROR << "Z'->ee reconstruction on AOD failed" << endreq;
    return StatusCode::SUCCESS;
  }
  
  return sc;

}

//////////////////////////////////////////////////////////////////////////////////
/// zprime on aod: called by execute()

StatusCode MCProduction::mcproduction() {

  MsgStream mLog( messageService(), name() );

  mLog << MSG::DEBUG << "mcproduction()" << endreq;

  StatusCode sc = StatusCode::SUCCESS;

  /// read the AOD electron container from persistecy storage
  const ElectronContainer* elecTES = 0;
  sc=m_storeGate->retrieve( elecTES, m_electronContainerName);
  if( sc.isFailure()  ||  !elecTES ) {
     mLog << MSG::WARNING
          << "No AOD electron container found in TDS"
          << endreq; 
     return sc;
  }  
  mLog << MSG::DEBUG << "ElectronContainer successfully retrieved" << endreq;

  const EventInfo *evt = 0;

  sc=m_storeGate->retrieve(evt, m_eventInfoKey);
  if( sc.isFailure()  || !evt ) {
    mLog << MSG::ERROR << "Error retrieving " << m_eventInfoKey << "!"<<endreq;
    return StatusCode::FAILURE;
  }
    
  int Event_Run_Number = evt->event_ID()->run_number();
  int Event_Event_Number = evt->event_ID()->event_number();
  int Event_Lumi_Block = evt->event_ID()->lumi_block();

  double mee = -1.0;
  int i = 0;

  // NEvents
  h_CutFlow->Fill(0);

  // NGRL
  h_CutFlow->Fill(1);

  if (1 == 1){

    // NVertices
    h_CutFlow->Fill(2);

    if (1 == 1){
  
      // NTrigger
      h_CutFlow->Fill(3);

    /// iterators over the container 
    ElectronContainer::const_iterator elecItr  = elecTES->begin();
    ElectronContainer::const_iterator elecItrE = elecTES->end();
   
    int ElecContainerSize = elecTES->size();

    if (ElecContainerSize >= 2){

    std::vector<ElectronContainer::const_iterator> GoodElectrons; 
    GoodElectrons.resize(ElecContainerSize);

    // Look for the most energetic electrons and sort them in decreasing energy
    for (; elecItr != elecItrE; ++elecItr) {
 
        GoodElectrons[i] = elecItr;
        i++;      
        
    }
 
    if (GoodElectrons.size() >= 2){

    //isEM Author Cut

    for (unsigned int a=0 ; a<GoodElectrons.size() ; a++){
    
      if ( (((*(GoodElectrons[a]))->author()) != 0) && (((*(GoodElectrons[a]))->author()) != 3) ){
    
        GoodElectrons.erase(GoodElectrons.begin()+a);
        a--;
      }
    }
  
    if ( GoodElectrons.size() >= 2 ){
  
      // NAuthor
      h_CutFlow->Fill(4);

      // Cluster Eta Cut

      for (unsigned int a=0 ; a<GoodElectrons.size() ; a++){

        double cleta = (*(GoodElectrons[a]))->eta();

        if ( !((fabs(cleta) < m_etaElecCut) && ((fabs(cleta) > 1.52) || fabs(cleta) < 1.37)) ){

          GoodElectrons.erase(GoodElectrons.begin()+a);
          a--;
	}
      }
      
      if ( GoodElectrons.size() >= 2 ){

        // NEta
        h_CutFlow->Fill(5);

        // Cluster Pt Cut

        for (unsigned int a=0 ; a<GoodElectrons.size() ; a++){

          double clpt = (*(GoodElectrons[a]))->pt();

          if ( !(clpt > m_ptElecCut) ){

            GoodElectrons.erase(GoodElectrons.begin()+a);
            a--;
	  }
	}
    
        if ( GoodElectrons.size() >= 2 ){

          // NPt
          h_CutFlow->Fill(6);

          // OTX Cut

          // MC or Data?
          if (m_MCorData == "MC"){
            Event_Run_Number = m_LatestOQ;
	  }

          for (unsigned int a=0 ; a<GoodElectrons.size() ; a++){

	    double cleta = (*(GoodElectrons[a]))->eta();
            double clphi = (*(GoodElectrons[a]))->phi();
            //int OTX = egammaOQ::checkOQCluster(Event_Run_Number,cleta,clphi,1);
            int OTX = 1;

            if ( OTX == 3 ){

              GoodElectrons.erase(GoodElectrons.begin()+a);
              a--;
	    }
	  }

          if ( GoodElectrons.size() >= 2 ){  
	 
            // NOTX
            h_CutFlow->Fill(7);

            // isEM Medium Cut

            for (unsigned int a=0 ; a<GoodElectrons.size() ; a++){

	      //bool isEMLoose = ((*(GoodElectrons[a]))->isem(egammaPID::Loose)==0);
	      bool isEMMedium = ((*(GoodElectrons[a]))->isem(egammaPID::ElectronMedium_WithTrackMatch)==0);
   	      //bool isEMTight = ((*(GoodElectrons[a]))->isem(egammaPID::Tight)==0);

              if ( !(isEMMedium) ){

                GoodElectrons.erase(GoodElectrons.begin()+a);
                a--;
	      }
	    }

            if ( GoodElectrons.size() >= 2 ){ 

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

            float TestPt1 = 0;
            float TestPt2 = 0;
            int firstrun = 0;

            std::vector<ElectronContainer::const_iterator> PairElectrons; 
            PairElectrons.resize(GoodElectrons.size());

            // Take Top Two Electrons to Make Pair
	    for (unsigned int a=0 ; a<GoodElectrons.size() ; a++ ) {

              float CandidatePt = ((*(GoodElectrons[a]))->pt());

              if ( CandidatePt > TestPt1 ){
                
                if ( firstrun == 1 ){
                  float Pair0pt = ((*(PairElectrons[0]))->pt());
                   
                  if ( Pair0pt > TestPt2 ){
                    PairElectrons[1] = PairElectrons[0];
                    TestPt2 = Pair0pt;
	  	  }

		}

                if ( firstrun == 0 ){
                  firstrun = 1;
		}

                PairElectrons[0] = GoodElectrons[a];
                TestPt1 = CandidatePt;

	      }
	      else if ( CandidatePt > TestPt2 ){
                PairElectrons[1] = GoodElectrons[a];
                TestPt2 = CandidatePt;
	      }
        
	    }

	    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
	    /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
            
         
            // Check Charge, Take 2 Top Electrons

            //if ( ((*(PairElectrons[0]))->charge() == -((*(PairElectrons[1]))->charge())) ){
	    if ( 0 == 0 ){
   
              h_CutFlow->Fill(11);

              // Normal Stuff

              TVector3 p1(1);
              p1.SetMag((*(PairElectrons[0]))->e());
	      //p1.SetEta((*(PairElectrons[0]))->eta());
              //p1.SetTheta((*(PairElectrons[0]))->theta());
              //p1.SetPhi((*(PairElectrons[0]))->phi());
              p1.SetPtEtaPhi((*(PairElectrons[0]))->pt(),(*(PairElectrons[0]))->eta(),(*(PairElectrons[0]))->phi());

	      TVector3 p2(1);
              p2.SetMag((*(PairElectrons[1]))->e());
              //p2.SetEta((*(PairElectrons[1]))->eta());
	      //p2.SetPhi((*(PairElectrons[1]))->phi());
              p2.SetPtEtaPhi((*(PairElectrons[1]))->pt(),(*(PairElectrons[1]))->eta(),(*(PairElectrons[1]))->phi());


              mee = sqrt(2*(p1.Mag() * p2.Mag()- p1*p2 ));


              float etafirst = (*(PairElectrons[0]))->eta();

              float Etfirst = ((*(PairElectrons[0]))->e())/cosh(etafirst);
              float Ptfirst = (*(PairElectrons[0]))->pt();
              float Etafirst = (*(PairElectrons[0]))->eta();

              float etasecond = (*(PairElectrons[1]))->eta();

              float Etsecond = ((*(PairElectrons[1]))->e())/cosh(etasecond);
              float Ptsecond = (*(PairElectrons[1]))->pt();
              float Etasecond = (*(PairElectrons[1]))->eta();
                 
              // Fill Pre Invariant Mass and Kinematics
              h_Pair_Pt->Fill(Ptfirst/1000, m_Scale);
              h_Pair_Et->Fill(Etfirst/1000, m_Scale);
              h_Pair_Eta->Fill(etafirst, m_Scale);

              h_Pair_Pt->Fill(Ptsecond/1000, m_Scale);
              h_Pair_Et->Fill(Etsecond/1000, m_Scale);
              h_Pair_Eta->Fill(etasecond, m_Scale);

              h_Pair_InvMass->Fill(mee/1000, m_Scale);

	    }
	    }
	  }
	}
      }
    }
    }
    }
    }
  }
  return StatusCode::SUCCESS;
}