#define CollectionTree_cxx
#include "CollectionTree.h"
#include <TH2.h>
#include <TStyle.h>
#include <stdio.h>
#include <TCanvas.h>
#include <vector>
#include <algorithm>
#include <Riostream.h>
#include <math.h>
#include <fstream.h>
#include "TCanvas.h"
#include "TH1.h"
#include "TPad.h"
#include "TFile.h"
#include "TStopwatch.h"
#include <string.h>


using namespace std;

class Jet {
public:
  Jet();
  
  Double_t JetWeight;
  Int_t    Label;
  Int_t PdgCode;
  Int_t Goodness;
  Float_t Eta;
  Float_t Phi;
  float pT;

  bool operator<(const Jet& jet) const { 
	return this->JetWeight>jet.JetWeight;

	}
};


Jet::Jet(){

  JetWeight=0;
  Label=0;
  PdgCode=0;
  Goodness=0;
  Eta=0;
  Phi=0;
  pT=0;


}

double  Calculate_DR(Jet a, Jet b){
  double Dr= pow ((a.Eta-b.Eta )*(a.Eta-b.Eta )+(a.Phi-b.Phi)*(a.Phi-b.Phi), 0.5);
  return Dr;
}




void PaintBin (TH1 *h, Int_t bin, Int_t color) 
{
   printf("%d %d %d\n", bin, color, h->GetBinContent(bin));
   TBox *b = new TBox(h->GetBinLowEdge(bin),
                      h->GetMinimum(),
                      h->GetBinWidth(bin)+h->GetBinLowEdge(bin),
                      h->GetBinContent(bin));
   b->SetFillColor(color);
   b->Draw();

}


void CollectionTree::Loop()
{

  FILE *Output_L;
  Double_t TotalShareDR =0 , TotalSharepT =0 ; 
  Double_t TotalShare =0; 
  Int_t  NtotJet=0;
  Int_t Ntotb=0;
  float DRref; 
  Double_t FractpTref;
  double DR;
  Int_t Switch1;
  float pTMax =0.0;

  


TStopwatch *st=new TStopwatch();
st->Start();


  std::cout << " ------------------------------- " << std::endl;
  std::cout << " Improving b-tagging efficiency:" << std::endl;
  std::cout << " ------------------------------- " << std::endl;

  int NumCols=1;
  std::vector<int> *a = new std::vector<int>[NumCols];
  std::vector<double> DRcollection; 
  std::vector<float> pTcollection; 
  std::vector<float> pTcollcleaned;
  std::vector<double> DRcollcleaned ;
  std::vector<float> ppT;
 
   std::cout << "Insert minimum angular separation to distinguish different jet: " << endl;
   std::cin >> DRref;
   if (DRref >= 8.0) {
   cerr <<"Too high value."<< endl;
   std::cout << "Insert minimum angular separation to distinguish different jet: " << endl;
   std::cin >> DRref;
 }

   Output_L = fopen( "Output_L_C.txt", "a" );

   if (Output_L == NULL) {
	std::cout << "gioWARNING: No output file has been opened." << std::endl;
        exit(1);
     }
   fprintf( Output_L, "\n·.·´¯`·.··.·´¯`·.··.·´¯`·.··.·´¯`·.··.·´¯`·.··.·´¯`·.·\n" );
   fprintf( Output_L, "Minimum angular separation to distinguish different jet: %f", DRref , "\n" ); 
 
  Int_t NbinsDR = 120;
  Float_t MaxHistoValueDR = 8.0;
  Int_t highlightedBinDR = static_cast<int>(DRref*NbinsDR/MaxHistoValueDR);

std::cout << "    Trasverse momentum" << std::endl;
  std::cout << "Insert N, where N is pTMax/N= pTMin to be considered as a jet: " << endl;
  std::cin >> FractpTref;


TH1F *h1 = new TH1F("histoDR","Angular Separation", NbinsDR , 0.0 , MaxHistoValueDR );

h1->GetXaxis()->SetTitle("DR");
h1->GetYaxis()->SetTitle("Counts");
h1->SetFillColor(38);
 Int_t NbinspT = 120;
  Float_t MaxHistoValuepT = 700000;


TH1F *h4 = new TH1F("histopT2","Transverse Momentum - Cut", NbinspT , 0.0 , MaxHistoValuepT );

h4->GetXaxis()->SetTitle("pT [MeV]");
h4->GetYaxis()->SetTitle("Counts");
h4->SetFillColor(20);


  if (fChain == 0) return;

  Long64_t nentries = fChain->GetEntriesFast();

  Long64_t nbytes = 0, nb = 0;

  int NentDR=0;
  int NentpT=0;
  int Nent=0;  
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    Long64_t ientry = LoadTree(jentry);
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;

    int NtotbJet=0;
    vector <Jet> jet;
    for (int i=0; i<PJet_N; i++) {
 
      Jet tmp;
      tmp.JetWeight= (*PJet_weight)[i];
      tmp.PdgCode=(*PJet_pdgId)[i];
      tmp.Label= i;
      tmp.Eta=(*PJet_eta)[i];
      tmp.Phi=(*PJet_phi)[i];
      tmp.pT=(*PJet_p_T)[i];
     
      jet.push_back(tmp);
    }


    for (int i=0;i<jet.size();i++){
      if (jet[i].PdgCode == 5) jet[i].Goodness = 1;
      else jet[i].Goodness = 0;
    
    }

    sort(jet.begin(),jet.end());

   

    for (int i=0;i<jet.size();i++){

      NtotbJet =  NtotbJet + jet[i].Goodness;
    }
      

    NtotJet = NtotJet + jet.size(); 


    int SumFourGoodness =0;
    int SumFourGoodness4=0;
 
    for (int i=0 ; i<jet.size(); i++){


    SumFourGoodness = SumFourGoodness + jet[i].Goodness;
    if (i<4) 
	SumFourGoodness4 = SumFourGoodness4 + jet[i].Goodness;     
    }
  
    for (int cols=0 ; cols< NumCols ; cols++){

    	a[cols].clear();
    }
    
         NumCols=1;

    a[0].push_back(jet[0].Label);

    for (int i=0; i<jet.size(); i++){
      if (i>0)
	for (int cols=0; cols<NumCols; cols++){
	  int b = a[cols].front();
	  DR = Calculate_DR(jet[i],jet[b]);	 

	  if (DR> DRref){
	    NumCols++;
	    
	    std::vector<int> *temp= new std::vector<int>[NumCols]; 
	    for (int k = 0 ; k < NumCols -1; k++) temp[k]=a[k]; 
	    a = temp;  
	    a[NumCols-1].push_back(jet[i].Label);
	      break;
	  }
	
	  if (DR <= DRref){
	    a[cols].push_back(jet[i].Label);
	    break;
	  }
	}
    }


  for (int i=0 ; i<jet.size(); i++){
	  for (int ii=0; ii<i ; ii++){
		DR = Calculate_DR(jet[i],jet[ii]);
		DRcollection.push_back(DR);
		}
	}


for (i=0; i<NumCols; i++){
	  for (int ii=0; ii<i ; ii++){
		DR = Calculate_DR(jet[a[i].front()],jet[a[ii].front()]);
		DRcollcleaned.push_back(DR);

		}

	}



Int_t SumGoodnessDR = 0;
Int_t SumFirst4GoodnessDR = 0;

for (i=0; i<NumCols; i++){

    SumGoodnessDR = SumGoodnessDR + jet[i].Goodness;
    if (i<4){
      SumFirst4GoodnessDR = SumFirst4GoodnessDR + jet[i].Goodness;
      }
}

     double ShareDR=0;
    if (SumFirst4GoodnessDR==4) ShareDR=1;
    if (SumGoodnessDR>=4) NentDR++; 
    TotalShareDR  = TotalShareDR+ShareDR ;


    double Share=0;
    if (SumFourGoodness4==4) Share=1;
    if (SumFourGoodness>=4) Nent++; 
   
    TotalShare  = TotalShare+Share ;
   
    Ntotb= Ntotb + NtotbJet;

for (int i=0; i<jet.size(); i++){
       pTcollection.push_back(jet[i].pT);
}

for (int i=0; i<jet.size(); i++){
	if (jet[i].pT > pTMax){
		pTMax=jet[i].pT;
	}
		
	else pTMax = pTMax;
}

 ppT.clear();

 Double_t Fract = pTMax/FractpTref;

       for (int i=0; i<jet.size(); i++){
          
          if (jet[i].pT > Fract) {
          pTcollcleaned.push_back(jet[i].pT);
          ppT.push_back(jet[i].Label);
	}
}


int SumGoodnesspT =0;
  int Sum4GoodnesspT=0;

   
for (int i=0 ; i<ppT.size(); i++){
    SumGoodnesspT =  SumGoodnesspT +jet[ppT[i]].Goodness;
  if (i <4) Sum4GoodnesspT = Sum4GoodnesspT + jet[ppT[i]].Goodness;
}	



double SharepT=0;
    
    if (Sum4GoodnesspT==4) SharepT=1;
    if (SumGoodnesspT>=4) NentpT++; 
   
    TotalSharepT  = TotalSharepT + SharepT ;    



    }//______________________________________-


TH1F *h3 = new TH1F("histoDR","Angular Separation - Cut ", NbinsDR , 0.0 , MaxHistoValueDR );

h3->GetXaxis()->SetTitle("#DeltaR");
h3->GetYaxis()->SetTitle("Counts");
h3->SetFillColor(38);



TH1F *h2 = new TH1F("histopT","Transverse Momentum", NbinspT , 0.0 , MaxHistoValuepT );

h2->GetXaxis()->SetTitle("pT [MeV]");
h2->GetYaxis()->SetTitle("Counts");
h2->SetFillColor(20);




for (int g=0; g<DRcollection.size(); g++){
	h1->Fill( DRcollection[g]);
}

for (int g=0; g<pTcollection.size(); g++){
	h2->Fill( pTcollection[g]);
}

for (int g=0; g<DRcollcleaned.size(); g++){
	h3->Fill( DRcollcleaned[g]);
}

  Double_t PercJetSeparatSelected =  (TotalShareDR*100/NentDR);
  Double_t PercJetpTSelected =  (TotalSharepT*100/NentpT);


st->Stop();
st->RealTime();
st->CpuTime();
st->Print(" "); 

for (int g=0; g<pTcollcleaned.size(); g++){ 
	h4->Fill( pTcollcleaned[g]);
}

TCanvas *Histogram1 = new TCanvas("Histogram1","Histogram1",1500,1100);

Histogram1->Divide(2,2);
Histogram1->cd(1);

h1->Draw();

for (int iBin=0; iBin<= highlightedBinDR; iBin++){
PaintBin (h1, iBin , kRed);
}

Histogram1->cd(2);
h2->Draw();
 
Int_t highlightedBinpT = static_cast<int>((pTMax/FractpTref)*NbinspT/MaxHistoValuepT);
for (int iBin=0; iBin<= highlightedBinpT; iBin++){
PaintBin (h2, iBin , kRed);
}

Histogram1->cd(3);
h3->Draw();

Histogram1->cd(4);
h4->Draw();


Histogram1->Print("Histo1.ps");


  std::cout << "?·.·´¯`·.·??·.·´¯`·.·?·.·´¯`·.·??·.·´¯`·.·??·.·´¯`·.·?·.·´¯`·.·??" << std::endl;
  std::cout << " ------------------------------- " << std::endl;
  std::cout << "        S U M M A R Y:  " << std::endl;
  std::cout << " ------------------------------- " << std::endl;
  std::cout << " Number of Trees in the Chain: " << nentries << std::endl;
  std::cout << " Number of event with >=4 b-jet:  " << Nent << std::endl;
  std::cout << " Total Number of jets: " << NtotJet << std::endl;
  std::cout << " Total number of b-jet (referring to PJet_pdgId) : " <<   Ntotb << std::endl;

  std::cout << " ------------------------------- " << std::endl;
  std::cout << "      E F F I C I E N C Y:  " << std::endl;
  std::cout << " ------------------------------- " << std::endl;
  std::cout << "\n ";
  std::cout << "1) Average Share (4 higher PJet_weight values = 4 b-jet) over all the Trees: " << (TotalShare*100/Nent) << "%" << std::endl;
  std::cout << "\n ";
  std::cout << "2) Average Share With Restriction on Direction: " <<   PercJetSeparatSelected << "%" << std::endl;
  std::cout << "\n ";
  std::cout << "3) Average Share With Restriction on Trasverse impulse: " <<  PercJetpTSelected << "%" << std::endl;

  fprintf( Output_L, "Average Share (4 higher PJet_weight values = 4 b-jet) over all the Trees: %f %%", (TotalShare*100/Nent) , "\n" );
  fprintf( Output_L, "Average Share With Restriction on Direction: %f %%", PercJetSeparatSelected , "\n" );
  fprintf( Output_L, "Average Share With Restriction on Trasverse impulse: %f %%", PercJetpTSelected , "\n" );
  fclose (Output_L);


}