#include <TVectorT.h>
#include <TROOT.h>
#include <TString.h>
#include <TFile.h>
#include <TGraph.h>
#include <TStyle.h>
#include <TH1D.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TLegend.h>
#include <TList.h>
#include <iostream>
#include <cstdio>   // needed for io
#include <vector>
#include <TSystem.h>

#include "/usr/local/include/fastjet/PseudoJet.hh"    
//#include "fastjet/PseudoJet.hh"

//using namespace fastjet;
//using namespace std;

#ifdef __CINT__

#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ nestedclasses;
#pragma link C++ class fastjet::PseudoJet+;
#pragma link C++ class std::vector<fastjet::PseudoJet>+;

#endif

void NewScaling(TH1D *h, float norma );

void dsigdpt(){

  gSystem->Load("/usr/local/lib/libfastjet.so");
  gSystem->Load("/usr/local/lib/libfastjettools.so");
  gSystem->AddIncludePath("/usr/local/include/fastjet");
 
  cout<<"entering "<<endl;  
  TFile *fdCT10 = new TFile("tryforphonly/ph_jet/ptm01/isol5GeV/ggdtest_n2CMSCT10nlo.root"); 
  TTree *t4 = (TTree*)fdCT10->Get("t2");

  TFile *foCT10 = new TFile("tryforphonly/ph_jet/ptm01/isol5GeV/ggotest_n2CMSCT10nlo.root"); 
  TTree *t5 = (TTree*)foCT10->Get("t2");

  
  TH1D *p86682= new TH1D("p86682","dsigma_lo/y_gamma",100,20,1000);

  //vector<fastjet::PseudoJet> input_particles;
  //
  
  Int_t iprov2,ntrack2,iprov3,ntrack3;
  
  Double_t e2[3],px2[3],py2[3],pz2[3];
  Double_t e3[3],px3[3],py3[3],pz3[3];
 
  
  Double_t ptj2[2],yj2[2];
  Double_t ptj3[2],yj3[2];
    

  Double_t pt2[3],y2[3],phi2[3];
  Double_t pt3[3],y3[3],phi3[3];
 
  Float_t pdf_weight2[1000];
  Float_t pdf_weight3[1000];
  
  Float_t weight2,weight3;

  // we get the value stored into the header for the normalisation

  TList* list2 = t4->GetUserInfo();
  TList* list3 = t5->GetUserInfo();
  
    
  TVectorT<float> &v2 = *(list2->At(0));
  TVectorT<float> &v3 = *(list3->At(0)); 

  
  float& nb_evt2 = v2[0];
  float& nb_evt3 = v3[0];
    
  float& xsec2 = v2[1];
  float& xsec3 = v3[1];
    

  float& sqrt_s2 = v2[2];
  float& sqrt_s3 = v3[2];

    
 
  float norma2 = xsec2/nb_evt2;
  float norma3 = xsec3/nb_evt3;
  


  t4->SetBranchAddress("iprov",&iprov2);
  t4->SetBranchAddress("ntrack",&ntrack2);
  t4->SetBranchAddress("energy",e2);
  t4->SetBranchAddress("px",px2);
  t4->SetBranchAddress("py",py2);
  t4->SetBranchAddress("pz",pz2);
  t4->SetBranchAddress("pdf_weight",pdf_weight2);


  t5->SetBranchAddress("iprov",&iprov3);
  t5->SetBranchAddress("ntrack",&ntrack3);
  t5->SetBranchAddress("energy",e3);
  t5->SetBranchAddress("px",px3);
  t5->SetBranchAddress("py",py3);
  t5->SetBranchAddress("pz",pz3);
  t5->SetBranchAddress("pdf_weight",pdf_weight3);


 
  
  //
  Int_t bin_pt = 100;
  Int_t bin_y = 100;
  Double_t pt_min = 40.;
  Double_t pt_max = 300.;
  Double_t y_min= -0.9;
  Double_t y_max= 0.9;
  Double_t bin_size_pt;

  bin_size_pt = (pt_max-pt_min)/(Double_t) bin_pt;
  //
  TH1D *hpt = new TH1D("pt","essai",bin_pt,pt_min,pt_max);
  //
  TH1D *hpt1 = new TH1D("pt1","essai",bin_pt,pt_min,pt_max);
  int nbins=8;
  double xval[9] = {40.0, 45.0, 50.0, 60.0, 70.0, 85.0, 100.0, 145.0, 300.0};

  TH1D *CT10hdcentral = new TH1D("CT10hdcentral","",nbins,xval);
  TH1D *CT10hocentral = new TH1D("CT10hocentral","",nbins,xval);
  TH1D *CT10hdmin = new TH1D("CT10hdmin","",nbins,xval);
  TH1D *CT10hdmax = new TH1D("CT10hdmax","",nbins,xval);
  TH1D *CT10homin = new TH1D("CT10homin","",nbins,xval);
  TH1D *CT10homax = new TH1D("CT10homax","",nbins,xval);
  
  
  CT10hdcentral->Sumw2();
  CT10hocentral->Sumw2();
  CT10hdmin->Sumw2();
  CT10hdmax->Sumw2();
  CT10homin->Sumw2();
  CT10homax->Sumw2();

  
    
  Int_t entries2 = (Int_t)t4->GetEntries();
  Int_t entries3 = (Int_t)t5->GetEntries();    
 
  Double_t pi,r12;
  Double_t r_kt;
  pi =3.141592653589793; 
  r_kt = 0.5;
  Bool_t merged;
  
  Double_t rpj1,rpj2,rpj;
  
  double max2, min2;
  for (Int_t i=0;i<entries2;i++) {
    t4->GetEntry(i);
    // Pt and y are built, 0 is the photon, 1 is the hard parton and 2 is the soft one
    for (Int_t j=0;j<ntrack2;j++) {     
      pt2[j] = sqrt(px2[j]*px2[j]+py2[j]*py2[j]);
      y2[j] = log( (e2[j]+pz2[j])/(e2[j]-pz2[j]) )*0.5;
      phi2[j] = atan2(py2[j],px2[j]);
      if (phi2[j] < 0.) {
        phi2[j] = phi2[j] + 2.*pi;
      }
      //  if(j>0) { input_particles.push_back(fastjet::PseudoJet(px2[j],py2[j],pz2[j],e2[j]) );}
    }
    /*  fastjet::JetAlgorithm jet_alg = fastjet::antikt_algorithm;
    double R = 0.5;
    fastjet::Strategy strategy = fastjet::Best;
    fastjet::JetDefinition jet_def(jet_alg,R,strategy);
    fastjet::ClusterSequence clust_seq(input_particles, jet_def);
    double ptmin=40.0;
    vector<fastjet::PseudoJet> inclusive_jets = sorted_by_pt(clust_seq.inclusive_jets(ptmin));
    cout << "Ran " << jet_def.description() << endl;
    printf("%5s %15s %15s %15s\n","jet #", "rapidity", "phi", "pt");

    for (unsigned int i = 0; i < inclusive_jets.size(); i++) {
      printf("%5u %15.8f %15.8f %15.8f\n",
	     i, inclusive_jets[i].rap(), inclusive_jets[i].phi(),
	     inclusive_jets[i].perp());
    }
    */
    weight2 = pdf_weight2[0];
    max2 = pdf_weight2[0];
    min2 = pdf_weight2[0];
    
    if (ntrack2 == 3) {
      r12 = sqrt((y2[1]-y2[2])*(y2[1]-y2[2])+(phi2[1]-phi2[2])*(phi2[1]-phi2[2]));
      if (r12 <= r_kt) {
	e2[1] = e2[1]+e2[2];
	px2[1] = px2[1]+px2[2];
	py2[1] = py2[1]+py2[2];
	pz2[1] = pz2[1]+pz2[2];
	phi2[1]=(phi2[1]+phi2[2])/2;
	e2[2] = 0.;
	px2[2] = 0.;
	py2[2] = 0.;
	pz2[2] = 0.;
	phi2[2]=0;
	
      }
    }
 
    for (Int_t j=1;j<=2;j++) {
      ptj2[j-1] =sqrt(px2[j]*px2[j]+py2[j]*py2[j]);
      if (e2[j] == 0.) {
        yj2[j-1] = 0.;
      }
      else {
        yj2[j-1] = log( (e2[j]+pz2[j])/(e2[j]-pz2[j]) )*0.5;
      }
    }
    
    for(int k=1;k<54;k++){
      if(pdf_weight2[k]>max2 && pdf_weight2[k]!=0){ max2 = pdf_weight2[k]; }
      if(pdf_weight2[k]<min2 && pdf_weight2[k]!=0){ min2 = pdf_weight2[k]; }
    }
    // fill histograms
   
      if(ptj2[0]>30 && ptj2[1]>30 && fabs(yj2[0])<1.5 ){ 
	CT10hdmax->Fill(pt2[0],max2);
	CT10hdmin->Fill(pt2[0],min2);
	CT10hdcentral->Fill(pt2[0],weight2);
      }
  }
  	    
  double max3, min3;
  for (Int_t i=0;i<entries3;i++) {
    t5->GetEntry(i);
    for(int l=0;l<ntrack3;l++){
      pt3[l] = sqrt(px3[l]*px3[l]+py3[l]*py3[l]);
      y3[l] = log( (e3[l]+pz3[l])/(e3[l]-pz3[l]) )*0.5;
      phi3[l] = atan2(py3[l],px3[l]);
      if (phi3[l] < 0.) {
	phi3[l] = phi3[l] + 2.*pi;
      }
    } 
    weight3 = pdf_weight3[0];
    max3 = pdf_weight3[0];
    min3 =   pdf_weight3[0];
    merged = false;
    
    if (ntrack3 == 3) {
      r12 = sqrt((y3[1]-y3[2])*(y3[1]-y3[2])+(phi3[1]-phi3[2])*(phi3[1]-phi3[2]));
      if (r12 <= r_kt) {
	merged=true;
	e3[1] = e3[1]+e3[2];
	px3[1] = px3[1]+px3[2];
	py3[1] = py3[1]+py3[2];
	pz3[1] = pz3[1]+pz3[2];
	phi3[1]= (phi3[1]+phi3[2])/2;
	e3[2] = 0.;
	px3[2] = 0.;
	py3[2] = 0.;
	pz3[2] = 0.;
	phi3[2]=0;
      }
    }  
    for (Int_t j=1;j<=2;j++) {
      ptj3[j-1] =sqrt(px3[j]*px3[j]+py3[j]*py3[j]);
      if (e3[j] == 0.) {
        yj3[j-1] = 0.;
      }
      else {
        yj3[j-1] = log( (e3[j]+pz3[j])/(e3[j]-pz3[j]) )*0.5;
      }
    }
    
    for(int k=1;k<54;k++){
      if(pdf_weight3[k]>max3 && pdf_weight3[k]!=0){ max3 = pdf_weight3[k]; }
      if(pdf_weight3[k]<min3 && pdf_weight3[k]!=0){ min3 = pdf_weight3[k]; }
    }
    
    // fill histograms
      if(ptj3[0]>30 && ptj3[1]>30 && fabs(yj3[0])<1.5 ){  	   

	CT10homax->Fill(pt3[0],max3);
	CT10homin->Fill(pt3[0],min3);
	CT10hocentral->Fill(pt3[0],weight3);
    }
  }
  
  cout<<"fourth step finished"<<endl;	  
    
  // normalisation of the histograms
   
   
    NewScaling(CT10hdmax,norma2);
    NewScaling(CT10hdmin,norma2);
    NewScaling(CT10hdcentral,norma2);

    NewScaling(CT10homax,norma3);
    NewScaling(CT10homin,norma3);
    NewScaling(CT10hocentral,norma3);
  
  cout<<"7th step finished"<<endl;	    
 
  CT10hdmax->Add(CT10homax);
  CT10hdmin->Add(CT10homin);
  CT10hdcentral->Add(CT10hocentral);
  
 
  p86682->SetLineColor(1);
  p86682->SetLineWidth(2);
  p86682->SetMarkerStyle(8);
  p86682->SetMarkerSize(1);
  cout<<"9th step finished"<<endl;	  
  
  //  TFile *fCMS = new TFile("rootFiles/ptCMSdata2.root");
  //  fCMS->ls();  
  //  c1->Draw();    
 
  CT10hdmin->SetLineColor(38);
  // CT10hdmin->Draw("same");
  CT10hdmax->SetLineColor(4);
  // CT10hdmax->Draw("same");
  CT10hdcentral->SetLineColor(36);
  // CT10hdcentral->Draw("same");
 CT10hdcentral->Draw();

  cout<<"11th step finished"<<endl;	  

  
  TLegend *leg0= new TLegend (0.2,0.75,0.6,0.88);

  
  // leg0->AddEntry(CT10hdmax,"CT10hdmax,NLO(Jetphox) CT10","pl");
  // leg0->AddEntry(CT10hdmin,"CT10hdmin,NLO(Jetphox) CT10","pl");
  leg0->AddEntry(CT10hdcentral,"CT10hdcentral,NLO(Jetphox) CT10","pl");
  leg0->AddEntry(p86682,"ph_jet Data 2011 #sqrt{s} =7 TeV","pl");
  leg0->SetFillColor(0);
  leg0->SetBorderSize(0);
  leg0->Draw();     
  cout<<"12th step finished"<<endl;	  
}
void NewScaling(TH1D *h, float norma){
  int nb_bins=8;
    Double_t xval[9] = {40.0, 45.0, 50.0, 60.0, 70.0, 85.0, 100.0, 145.0, 300.0};

  
  Double_t pt_bin_size[11]={0};
  
  for (int l=0;l<=(nb_bins-1);l++){
    pt_bin_size[l]=(xval[l+1]-xval[l]);
  }
  Double_t initial_bin_content,final_bin_content,initial_bin_err,final_bin_err;    
  for(int bin=1;bin<=nb_bins;bin++){
      
    initial_bin_content = 0;
    final_bin_content = 0;
    initial_bin_err=0;
    final_bin_err=0;
      
    initial_bin_content = h->GetBinContent(bin);
    initial_bin_err = h->GetBinError(bin);
      final_bin_content = ((initial_bin_content*norma)/(pt_bin_size[bin-1]));
      final_bin_err = ((initial_bin_err*norma)/(pt_bin_size[bin-1]));
      h->SetBinContent(bin,final_bin_content);
      h->SetBinError(bin,final_bin_err);
  }
}