// user include files
#include "FWCore/Framework/interface/Frameworkfwd.h"
#include "FWCore/Framework/interface/EDAnalyzer.h"

#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"

#include "DataFormats/Candidate/interface/Candidate.h"
#include "DataFormats/HcalRecHit/interface/HcalRecHitCollections.h"
#include "DataFormats/Common/interface/EDProduct.h"
#include "DataFormats/Common/interface/Ref.h"

#include "PhysicsTools/UtilAlgos/interface/TFileService.h"
#include "SimDataFormats/HepMCProduct/interface/HepMCProduct.h"

#include "HepMC/GenEvent.h"
#include "HepMC/HeavyIon.h"

#include <TMath.h>
#include <TH2D.h>
#include <TH1D.h>
#include <TFile.h>

using namespace std;

//
// class decleration
//

class HFAnalyzer : public edm::EDAnalyzer {
   public:
      explicit HFAnalyzer(const edm::ParameterSet&);
      ~HFAnalyzer();


   private:
      virtual void beginJob(const edm::EventSetup&) ;
      virtual void analyze(const edm::Event&, const edm::EventSetup&);
      virtual void endJob() ;

      // ----------member data ---------------------------

        double eta;
        double phi;
        double energy;
        double et;
        int ch;
        double eFwd;
        double eHF;
        double eHF1;
        double eHF2;
        double etHF1;
        double etHF2;
        double r;
        float b;
        double eCas;
        double rCas;
        double etCas;
        double E;
        double Et;


        TFile *outfile;

        TH2D *eHistE;
        TH2D *eHistHF;
        TH2D *eHistEt;
        TH2D *eHistF;
        TH2D *eHistR;
        TH2D *eCastE;
        TH2D *eCastEt;
        TH2D *eCastR;
        TH2D *eMidE;
        TH2D *eMidEt;
        TH2D *eMidR;
        TH1D *eHistb;
        TH1D *eHistDistb;
        TH1D *eHistDistb1;
        TH1D *eHistDistb2;
        TH1D *eHistDistb3;
        TH1D *eHisteta;
        TH1D *eHisteta1;
        TH1D *eHisteta2;
        TH2D *eHistK;
        TH2D *eHistL;
        TH1D *eHistphi;
        TH1D *eHistphi1;
        TH1D *eHistphi2;
        TH1D *eHistNch;
        TH1D *eHistetb;
        TH1D *eHistetb1;
        TH1D *eHistetb2;
        TH1D *eHistetb3;
        TH1D *eHistetb4;
        TH1D *eHistetb5;
        TH1D *eHistetb6;
        TH1D *eHistetb7;
        TH1D *eHisteHF;
        TH1D *eHistNchen;
        TH1D *eHistNchtrv;
        TH1D *eHistNchenphi;
        TH1D *eHistNchphi;
        TH1D *eHistNchtrvphi;

};

//
// constants, enums and typedefs
//

//
// static data member definitions
//

//
// constructors and destructor
//
HFAnalyzer::HFAnalyzer(const edm::ParameterSet& iConfig)
{
   //now do what ever initialization is needed

}


HFAnalyzer::~HFAnalyzer()
{

   // do anything here that needs to be done at desctruction time
   // (e.g. close files, deallocate resources etc.)

}


//
// member functions
//

// ------------ method called to for each event  ------------
void
HFAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup)
{
   using namespace edm;
   using namespace reco;
   using namespace HepMC;

   ////////////////// impact parameter //////////////////

   Handle<HepMCProduct> hepmc;

   iEvent.getByLabel("source",hepmc);

   const GenEvent *ev = hepmc->GetEvent();

   HeavyIon *hi;

   hi = ev->heavy_ion();

   b = hi->impact_parameter();

  
  /////// Energy Calculation - Generated Particles /////////

   eFwd = 0;
   eHF = 0;
   eHF1 = 0;
   etHF1 = 0;
   eHF2 = 0;
   etHF2 = 0;
   etFwd = 0;
   r = 0;
   eCas = 0;
   rCas = 0;
   etCas = 0;
   E = 0;
   Et = 0;
 

        Handle<CandidateCollection> genParticles;
        iEvent.getByLabel("genParticleCandidates", genParticles);
        for(size_t ipar = 0; ipar < genParticles->size(); ++ ipar) {
           const Candidate & p = (*genParticles)[ ipar ];


        eta = p.eta();
        phi = p.phi();
        energy = p.energy();
        et = p.et();
        ch = p.charge();
   
         
        if(p.numberOfDaughters() == 0){


 //////////// HF ////////////////////////

        if((eta*eta<25.) && (eta*eta>9.)){

           eFwd = eFwd + p.energy();
           etFwd = etFwd + p.et();
           r = eFwd / etFwd; }

////////////// CASTOR  /////////////////

        if((eta*eta<41.) && (eta*eta>27.)){

           eCas = eCas + p.energy();
           etCas = etCas + p.et();
           rCas = eCas / etCas; }

///////////// HF1 & HF2 //////////////////

        if(eta<5. && eta>3.){

           eHF1 = eHF1 + p.energy();
           etHF1 = etHF1 + p.et(); }


        if(eta<-3. && eta>-5.){

           eHF2 = eHF2 + p.energy();
           etHF2 = etHF2 + p.et(); }

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

        if(eta*eta<9. && eta*eta>1.){
        
            E = E + p.energy();
            Et = Et + p.et();}

 //////////// For Charged Particles //////////////////

        if(ch != 0) {
         eHistNch->Fill(eta);
         eHistNchen->Fill(eta,energy);
         eHistNchtrv->Fill(eta,et);
         eHistNchenphi->Fill(phi,energy);
         eHistNchphi->Fill(phi);
         eHistNchtrvphi->Fill(phi,et);

}

 ///////// dEt/deta for different impact parameters /////////////

        if(b>0. && b<=2.){
        eHistetb->Fill(eta,et);}

        if(b>2. && b<=4.) {
        eHistetb1->Fill(eta,et);}

        if(b>4. && b<=6.) {
        eHistetb2->Fill(eta,et);}

        if(b>6. && b<=8.) {
        eHistetb3->Fill(eta,et);}
 
        if(b>8. && b<=10.) {
        eHistetb4->Fill(eta,et);}

        if(b>10. && b<=12.) {
        eHistetb5->Fill(eta,et);}

        if(b>12. && b<=14.) {
        eHistetb6->Fill(eta,et);}

        if(b>14. && b<=16.) {
        eHistetb7->Fill(eta,et);}

}
  }

//////// Energy deposition in HF ////////////////

      Handle<HFRecHitCollection> hits;
           iEvent.getByLabel("hfreco",hits);
        for(size_t ihit = 0; ihit<hits->size(); ++ ihit){
           const HFRecHit & rechit = (*hits)[ ihit ];
        eHF = eHF + rechit.energy();
        }

        if(eHF<=103500. && eHF>=100000.){
        eHistDistb->Fill(b); }

        if(eHF<=73500. && eHF>=70000.){
        eHistDistb1->Fill(b); }

        if(eHF<=33500. && eHF>=30000.){
        eHistDistb2->Fill(b); }

        if(eHF<=3500. && eHF>=0.){
        eHistDistb3->Fill(b);}



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

   cout<<"GenEvent obtained"<<endl;

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

  //// Filling Histograms ////////////

    eHistE->Fill(b,eFwd);
    eHistHF->Fill(b,eHF);
    eHistEt->Fill(b,etFwd);
    eHistF->Fill(eFwd,eHF);
    eHistK->Fill(eHF1,eHF2);
    eHistL->Fill(etHF1,etHF2);
    eHistR->Fill(b,r);
    eHistb->Fill(b);
    eHisteHF->Fill(eHF);
    eCastE->Fill(b,eCas);
    eCastEt->Fill(b,etCas);
    eCastR->Fill(b,rCas);
    eMidE->Fill(b,E);
    eMidEt->Fill(b,Et);



#ifdef THIS_IS_AN_EVENTSETUP_EXAMPLE
   ESHandle<SetupData> pSetup;
   iSetup.get<SetupRecord>().get(pSetup);
#endif
    }


// ------------ method called once each job just before starting event loop  ------------
void
HFAnalyzer::beginJob(const edm::EventSetup&)
{

  outfile = new TFile("HFAnalyzer.root","RECREATE");
  eHistE = new TH2D("eHistE","Forward Particle Energy",300,0.,15.,1000,0.,300000.);
  eHistHF = new TH2D("eHistHF","HF Energy",100,0.,15.,200,0.,250000.);
  eHisteHF = new TH1D("eHisteHF","HF energy Dist",100,0.,250000.);
  eHistEt = new TH2D("eHistEt","Forward Particle Transverse Energy",300,0.,15.,1000,0.,15000.);
  eHistF = new TH2D("eHistF","Forward Particle Energy-HF",500,0.,350000.,500,0.,350000.);
  eHistR = new TH2D("eHistR","Ratio Between Fwd&etFwd",100,0.,15.,100,0.,50.);
  eHistK = new TH2D("eHistK","eHF1-eHF2",100,0.,250000.,100,0.,250000.);
  eHistL = new TH2D("eHistL","etHF1-etHF2",100,0.,10000.,100,0.,10000.);
  eHistb = new TH1D("eHistb","b Dist",100,0.,15.);
  eHistDistb = new TH1D("eHistDistb","b Dist 103.5>E>100",100,0.,15.);
  eHistDistb1 = new TH1D("eHistDistb1","b Dist 73.5>E>70",100,0.,15.);
  eHistDistb2 = new TH1D("eHistDistb2","b Dist 33.5>E>30",100,0.,15.);
  eHistDistb3 = new TH1D("eHistDistb3","b Dist 3.5>E>0",100,0.,15.);
  eHistNch = new TH1D("eHistNch","Charged Particle Dist",200,-10.,10.);
  eHistetb = new TH1D("eHistetb","dEt/d b=0-2",100,-10.,10.);
  eHistetb1 = new TH1D("eHistetb1","dEt/d b=2-4",100,-10.,10.);
  eHistetb2 = new TH1D("eHistetb2","dEt/d b=4-6",100,-10.,10.);
  eHistetb3 = new TH1D("eHistetb3","dEt/d b=6-8",100,-10.,10.);
  eHistetb4 = new TH1D("eHistetb4","dEt/d b=8-10",100,-10.,10.);
  eHistetb5 = new TH1D("eHistetb5","dEt/d b=10-12",100,-10.,10.);
  eHistetb6 = new TH1D("eHistetb6","dEt/d b=12-14",100,-10.,10.);
  eHistetb7 = new TH1D("eHistetb7","dEt/d b=14-16",100,-10.,10.);
  eCastE = new TH2D("eCastE","CASTOR Energy",300,0.,15.,1000,0.,500000.);
  eCastEt = new TH2D("eCastEt","CASTOR Transverse Energy",300,0.,15.,1000,0.,50000.);
  eCastR = new TH2D("eCastR","Ratio between E&Et",100,0.,15.,200,0.,500.);
  eMidE = new TH2D("eMidE","Mid Energy",100,0.,15.,1000,0.,200000.);
  eMidEt = new TH2D("eMidEt","Mid Trs Energy",100,0.,15.,1000,0.,25000.);
  eHistNchen = new TH1D("eHistNchen","Ch Energy eta",100,-10.,10.);
  eHistNchtrv = new TH1D("eHistNchtrv","Ch Trs En eta",100,-10.,10.);
  eHistNchenphi = new TH1D("eHistNchenphi","Ch Energy phi",100,-(TMath::Pi()),TMath::Pi());
  eHistNchphi = new TH1D("eHistNchphi","Ch part dist phi",100,-(TMath::Pi()),TMath::Pi());
  eHistNchtrvphi = new TH1D("eHistNchtrvphi","Ch trs En phi",100,-(TMath::Pi()),TMath::Pi());
   
}  

      // ------------ method called once each job just after ending the event loop  ------------
void
HFAnalyzer::endJob() {

outfile->Write();
outfile->Close();
return;
}

//define this as a plug-in
DEFINE_FWK_MODULE(HFAnalyzer);