#define UFO_jets_cxx
// root> T->Process("UFO_jets.C")
// root> T->Process("UFO_jets.C","some options")
// root> T->Process("UFO_jets.C+")
//


#include "UFO_jets.h"
#include <TH2.h>
#include <TStyle.h>
#include <iostream>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <TCanvas.h>


using namespace std;

Int_t fNumberOfEvents; /// Total number of events/
Double_t GeV = 0.001;
Int_t N_hist = 2;  //// Number of pT bins !!!!!!!

void UFO_jets::Begin(TTree * /*tree*/)
{
   TString option = GetOption();
   
   fNumberOfEvents = 0 ;
   cout<<" ************    INITIALIZING    *************"<<endl;

   cout<<option<<endl;
   
   ////vector of histos in pT bins
   h_CSSKPFlow_Trim_jet10_D2_bin = 0 ;
   h_CSSKUFO___Trim_jet10_D2_bin = 0 ;
   h_EMPFlow___Trim_jet10_D2_bin = 0 ;
   h_LCTopo____Trim_jet10_D2_bin = 0 ;
   h_LCUFO_____Trim_jet10_D2_bin = 0 ;
   h_TCC_______Trim_jet10_D2_bin = 0 ;
   h_Truth_____Trim_jet10_D2_bin = 0 ;
   
}

void UFO_jets::SlaveBegin(TTree * /*tree*/)
{
   TString option = GetOption();

   //***************** Histogram definitions *********************************************************

   ///// Vector of D2 histos 
   h_CSSKPFlow_Trim_jet10_D2_bin = new vector<TH1D*>;
   h_CSSKUFO___Trim_jet10_D2_bin = new vector<TH1D*>;
   h_EMPFlow___Trim_jet10_D2_bin = new vector<TH1D*>;
   h_LCTopo____Trim_jet10_D2_bin = new vector<TH1D*>;
   h_LCUFO_____Trim_jet10_D2_bin = new vector<TH1D*>;
   h_TCC_______Trim_jet10_D2_bin = new vector<TH1D*>;
   h_Truth_____Trim_jet10_D2_bin = new vector<TH1D*>;

      TString n_CSSKPFlow_Trim_jet10_D2_bin, n_CSSKUFO___Trim_jet10_D2_bin, n_EMPFlow___Trim_jet10_D2_bin, n_LCTopo____Trim_jet10_D2_bin, n_LCUFO_____Trim_jet10_D2_bin, n_TCC_______Trim_jet10_D2_bin, n_Truth_____Trim_jet10_D2_bin ;

   for (Int_t i=0; i < N_hist ; i++)
     {
       n_CSSKPFlow_Trim_jet10_D2_bin.Form("h_CSSKPFlow_Trim_jet10_D2_bin_%d",i);
       n_CSSKUFO___Trim_jet10_D2_bin.Form("h_CSSKUFO___Trim_jet10_D2_bin_%d",i);
       n_EMPFlow___Trim_jet10_D2_bin.Form("h_EMPFlow___Trim_jet10_D2_bin_%d",i);
       n_LCTopo____Trim_jet10_D2_bin.Form("h_LCTopo____Trim_jet10_D2_bin_%d",i);
       n_LCUFO_____Trim_jet10_D2_bin.Form("h_LCUFO_____Trim_jet10_D2_bin_%d",i);
       n_TCC_______Trim_jet10_D2_bin.Form("h_TCC_______Trim_jet10_D2_bin_%d",i);
       n_Truth_____Trim_jet10_D2_bin.Form("h_Truth_____Trim_jet10_D2_bin_%d",i);
    
       h_CSSKPFlow_Trim_jet10_D2_bin->push_back(new TH1D( n_CSSKPFlow_Trim_jet10_D2_bin.Data() ,"", 1000,0.0,100));
       h_CSSKUFO___Trim_jet10_D2_bin->push_back(new TH1D( n_CSSKUFO___Trim_jet10_D2_bin.Data() ,"", 1000,0.0,100));
       h_EMPFlow___Trim_jet10_D2_bin->push_back(new TH1D( n_EMPFlow___Trim_jet10_D2_bin.Data() ,"", 1000,0.0,100));
       h_LCTopo____Trim_jet10_D2_bin->push_back(new TH1D( n_LCTopo____Trim_jet10_D2_bin.Data() ,"", 1000,0.0,100));
       h_LCUFO_____Trim_jet10_D2_bin->push_back(new TH1D( n_LCUFO_____Trim_jet10_D2_bin.Data() ,"", 1000,0.0,100));
       h_TCC_______Trim_jet10_D2_bin->push_back(new TH1D( n_TCC_______Trim_jet10_D2_bin.Data() ,"", 1000,0.0,100));
       h_Truth_____Trim_jet10_D2_bin->push_back(new TH1D( n_Truth_____Trim_jet10_D2_bin.Data() ,"", 1000,0.0,100));

       fOutput->Add(h_CSSKPFlow_Trim_jet10_D2_bin->at(i));
       fOutput->Add(h_CSSKUFO___Trim_jet10_D2_bin->at(i));
       fOutput->Add(h_EMPFlow___Trim_jet10_D2_bin->at(i));
       fOutput->Add(h_LCTopo____Trim_jet10_D2_bin->at(i));
       fOutput->Add(h_LCUFO_____Trim_jet10_D2_bin->at(i));
       fOutput->Add(h_TCC_______Trim_jet10_D2_bin->at(i));
       fOutput->Add(h_Truth_____Trim_jet10_D2_bin->at(i));
     }
}

Bool_t UFO_jets::Process(Long64_t entry)
{
  fReader.SetLocalEntry(entry); // For a TChain

  //// Weight
  Double_t Lumi_0 = 1. ;
  Double_t XSecWeight = UFO_jets::GetXSxFiltEff(*mcChannelNumber) *  Lumi_0 ;
  Double_t WEIGHT = XSecWeight * EventInfoAuxDyn_mcEventWeights[0] ;

  //////////   Histogram filling ///////////////

  //// Histograms in pT bins
  Double_t Thuth_pT_min = 0., Thuth_pT_max = 2500. ;
  Double_t Thuth_pT_bin = ( Thuth_pT_max - Thuth_pT_min ) / N_hist ;
  Double_t Thuth_pT_cut = Thuth_pT_min + Thuth_pT_bin;
  for (Int_t i=0; i < N_hist ; i++)
    {
      for (int iCSSKPFlow_jet = 0, nCSSKPFlow_jet = (AntiKt10CSSKPFlowTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize()<2)? AntiKt10CSSKPFlowTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize():2 ; iCSSKPFlow_jet < nCSSKPFlow_jet; ++iCSSKPFlow_jet)
	{
	  if (    AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iCSSKPFlow_jet] * GeV >= Thuth_pT_min
	       && AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iCSSKPFlow_jet] * GeV <= Thuth_pT_cut)
	    {
	     h_CSSKPFlow_Trim_jet10_D2_bin -> at(i) -> Fill( AntiKt10CSSKPFlowTrimmedPtFrac5SmallR20JetsAuxDyn_D2[iCSSKPFlow_jet], WEIGHT );
	    }
	}
      for (int iCSSKUFO_jet = 0, nCSSKUFO_jet = (AntiKt10CSSKUFOTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize()<2)? AntiKt10CSSKUFOTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize():2 ; iCSSKUFO_jet < nCSSKUFO_jet; ++iCSSKUFO_jet)
	{
	  if (    AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iCSSKUFO_jet] * GeV >= Thuth_pT_min
	       && AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iCSSKUFO_jet] * GeV <= Thuth_pT_cut)
	    {
	     h_CSSKUFO___Trim_jet10_D2_bin -> at(i) -> Fill( AntiKt10CSSKUFOTrimmedPtFrac5SmallR20JetsAuxDyn_D2[iCSSKUFO_jet], WEIGHT );
	    }
	}
      for (int iEMPFlow_jet = 0, nEMPFlow_jet = (AntiKt10EMPFlowTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize()<2)? AntiKt10EMPFlowTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize():2 ; iEMPFlow_jet < nEMPFlow_jet; ++iEMPFlow_jet)
	{
	  if (    AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iEMPFlow_jet] * GeV >= Thuth_pT_min
	       && AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iEMPFlow_jet] * GeV <= Thuth_pT_cut)
	    {
	     h_EMPFlow___Trim_jet10_D2_bin -> at(i) -> Fill( AntiKt10EMPFlowTrimmedPtFrac5SmallR20JetsAuxDyn_D2[iEMPFlow_jet], WEIGHT );
	    }
	}
      for (int iLCTopo_jet = 0, nLCTopo_jet = (AntiKt10LCTopoTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize()<2)? AntiKt10LCTopoTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize():2 ; iLCTopo_jet < nLCTopo_jet; ++iLCTopo_jet)
	{
	  if (    AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iLCTopo_jet] * GeV >= Thuth_pT_min
	       && AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iLCTopo_jet] * GeV <= Thuth_pT_cut)
	    {
	     h_LCTopo____Trim_jet10_D2_bin -> at(i) -> Fill( AntiKt10LCTopoTrimmedPtFrac5SmallR20JetsAuxDyn_D2[iLCTopo_jet], WEIGHT );
	    }
	}
      for (int iLCUFO_jet = 0, nLCUFO_jet = (AntiKt10LCUFOTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize()<2)? AntiKt10LCUFOTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize():2 ; iLCUFO_jet < nLCUFO_jet; ++iLCUFO_jet)
	{
	  if (    AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iLCUFO_jet] * GeV >= Thuth_pT_min
	       && AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iLCUFO_jet] * GeV <= Thuth_pT_cut)
	    {
	     h_LCUFO_____Trim_jet10_D2_bin -> at(i) -> Fill( AntiKt10LCUFOTrimmedPtFrac5SmallR20JetsAuxDyn_D2[iLCUFO_jet], WEIGHT );
	    }
	}
      for (int iTCC_jet = 0, nTCC_jet = (AntiKt10TCCTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize()<2)? AntiKt10TCCTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize():2 ; iTCC_jet < nTCC_jet; ++iTCC_jet)
	{
	  if (    AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iTCC_jet] * GeV >= Thuth_pT_min
	       && AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iTCC_jet] * GeV <= Thuth_pT_cut)
	    {
	     h_TCC_______Trim_jet10_D2_bin -> at(i) -> Fill( AntiKt10TCCTrimmedPtFrac5SmallR20JetsAuxDyn_D2[iTCC_jet], WEIGHT );
	    }
	}
      for (int iTruth_jet = 0, nTruth_jet = (AntiKt10TruthTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize()<2)? AntiKt10TruthTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize():2 ; iTruth_jet < nTruth_jet; ++iTruth_jet)
	{
	  if (    AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iTruth_jet] * GeV >= Thuth_pT_min
	       && AntiKt10TruthTrimmedPtFrac5SmallR20JetsAux_pt[iTruth_jet] * GeV <= Thuth_pT_cut)
	    {
	     h_Truth_____Trim_jet10_D2_bin -> at(i) -> Fill( AntiKt10TruthTrimmedPtFrac5SmallR20JetsAuxDyn_D2[iTruth_jet], WEIGHT );
	    }
	}
      Thuth_pT_min += Thuth_pT_bin ;
      Thuth_pT_cut += Thuth_pT_bin;
    } // for loop on  N_hist

  /////////// PRINTING
  if(fNumberOfEvents%10000==0.)
    {
      cout<<" ------------------------------------------------------------"<<endl;
      cout<<" ************   Event #"<<fNumberOfEvents<< "    *************"<<endl;
      cout<<"mcChannelNumber: " << *mcChannelNumber <<endl;
      cout<<"MC weight: " << EventInfoAuxDyn_mcEventWeights[0] <<endl;
       cout<<"XSecWeight: " << XSecWeight <<endl;
      cout<<"WEIGHT: " << WEIGHT <<endl;
      cout<<" * JETS *"<<endl;
      cout<<"- CSSKPFlow  Jets: " << AntiKt10CSSKPFlowTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize() <<endl;
      cout<<"- CSSKUFO  Jets: " << AntiKt10CSSKUFOTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize() <<endl;
      cout<<"- EMPFlow  Jets: " <<  AntiKt10EMPFlowTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize() <<endl;
      cout<<"- LCTopo  Jets: " << AntiKt10LCTopoTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize() <<endl;
      cout<<"- LCUFO  Jets: " << AntiKt10LCUFOTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize() <<endl;
      cout<<"- TCC  Jets: " << AntiKt10TCCTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize() <<endl;
      cout<<"- Truth  Jets: " << AntiKt10TruthTrimmedPtFrac5SmallR20JetsAuxDyn_D2.GetSize() <<endl;
      cout<<" ------------------------------------------------------------"<<endl;
     }

  ++fNumberOfEvents;
  
   return kTRUE;
}

void UFO_jets::SlaveTerminate()
{
}

void UFO_jets::Terminate()
{
  cout<<"******** PER EVENT CALCULATION ENDED ******** "<<endl;
  cout<<"TOTAL NUMBER OF EVENTS (direct): "<<fNumberOfEvents<<endl;
  TString option = GetOption();

  /// ***************** Output file ***************
  TFile f1(option,"RECREATE");
  
  TString n_CSSKPFlow_Trim_jet10_D2_bin, n_CSSKUFO___Trim_jet10_D2_bin, n_EMPFlow___Trim_jet10_D2_bin, n_LCTopo____Trim_jet10_D2_bin, n_LCUFO_____Trim_jet10_D2_bin, n_TCC_______Trim_jet10_D2_bin, n_Truth_____Trim_jet10_D2_bin ;
   for (Int_t i=0; i < N_hist ; i++)
    {
      n_CSSKPFlow_Trim_jet10_D2_bin.Form("h_CSSKPFlow_Trim_jet10_D2_bin_%d",i);
      n_CSSKUFO___Trim_jet10_D2_bin.Form("h_CSSKUFO___Trim_jet10_D2_bin_%d",i);
      n_EMPFlow___Trim_jet10_D2_bin.Form("h_EMPFlow___Trim_jet10_D2_bin_%d",i);
      n_LCTopo____Trim_jet10_D2_bin.Form("h_LCTopo____Trim_jet10_D2_bin_%d",i);
      n_LCUFO_____Trim_jet10_D2_bin.Form("h_LCUFO_____Trim_jet10_D2_bin_%d",i);
      n_TCC_______Trim_jet10_D2_bin.Form("h_TCC_______Trim_jet10_D2_bin_%d",i);
      n_Truth_____Trim_jet10_D2_bin.Form("h_Truth_____Trim_jet10_D2_bin_%d",i);

      TH1D* outCSSKPFlow_Trim_jet10_D2_bin=dynamic_cast<TH1D*>(fOutput->FindObject(n_CSSKPFlow_Trim_jet10_D2_bin));
      TH1D* outCSSKUFO___Trim_jet10_D2_bin=dynamic_cast<TH1D*>(fOutput->FindObject(n_CSSKUFO___Trim_jet10_D2_bin));
      TH1D* outEMPFlow___Trim_jet10_D2_bin=dynamic_cast<TH1D*>(fOutput->FindObject(n_EMPFlow___Trim_jet10_D2_bin));
      TH1D* outLCTopo____Trim_jet10_D2_bin=dynamic_cast<TH1D*>(fOutput->FindObject(n_LCTopo____Trim_jet10_D2_bin));
      TH1D* outLCUFO_____Trim_jet10_D2_bin=dynamic_cast<TH1D*>(fOutput->FindObject(n_LCUFO_____Trim_jet10_D2_bin));
      TH1D* outTCC_______Trim_jet10_D2_bin=dynamic_cast<TH1D*>(fOutput->FindObject(n_TCC_______Trim_jet10_D2_bin));
      TH1D* outTruth_____Trim_jet10_D2_bin=dynamic_cast<TH1D*>(fOutput->FindObject(n_Truth_____Trim_jet10_D2_bin));

      outCSSKPFlow_Trim_jet10_D2_bin-> Write() ;
      outCSSKUFO___Trim_jet10_D2_bin-> Write() ;
      outEMPFlow___Trim_jet10_D2_bin-> Write() ;
      outLCTopo____Trim_jet10_D2_bin-> Write() ;
      outLCUFO_____Trim_jet10_D2_bin-> Write() ;
      outTCC_______Trim_jet10_D2_bin-> Write() ;
      outTruth_____Trim_jet10_D2_bin-> Write() ;
    }

  
  f1.Close();
  
  cout<<"Saved output in "<< option <<" yay! :D "<<endl;
}

double UFO_jets::GetXSxFiltEff(int mc_id ) //  xsection * gen_eff / Events // (Xsec[pb]  * kFact * filtereff)
{
  //jetjet
  if(mc_id==361022) return 1.;//2.4334E+06 * 3.3264E-04 / 30000 ;
  if(mc_id==361023) return 2.6454E+04 * 3.1953E-04 / 30000 ;
  if(mc_id==361024) return 2.5464E+02 * 5.3009E-04 / 30000 ;
  if(mc_id==361025) return 4.5536E+00 * 9.2325E-04 / 30000 ;
  if(mc_id==361026) return 2.5752E-01 * 9.4016E-04 / 30000 ;
  if(mc_id==361027) return 1.6214E-02 * 3.9282E-04 / 30000 ;
  if(mc_id==361028) return 6.2505E-04 * 1.0162E-02 / 30000 ;
  if(mc_id==361029) return 1.; //1.9640E-05 * 1.2054E-02 / 30000 ;
  if(mc_id==361030) return 1.1961E-06 * 5.8935E-03 / 30000 ;
  if(mc_id==361031) return 4.2260E-08 * 2.7015E-03 / 30000 ;
  if(mc_id==361032) return 1.0370E-09 * 4.2502E-04 / 30000 ;
  else return 1.;

}