#include <iostream>
#include <fstream>
#include <sstream>
#include <string>
#include <iomanip>
#include "TF1.h"
#include "TH1.h"
#include "TH2.h"
#include "TRandom3.h"
#include "TGraph.h"
#include "TGraphErrors.h"
#include "TFile.h"
#include "TVirtualFitter.h"
#include "MiniT.h"
#include "TCutG.h"
#include "TMatrixD.h"
using namespace std;
int Npar = 2;
//TFile ff("testWm_et20.root");
//TFile ff("../etmissb/data_noetcut.root");
//TFile ff("../etmissb/data_noetcut.root");
TFile ff("Wmunu_data.root");
//TFile ff("../etmiss2/dataWenu_noetcut_sumet.root");
//TFile gg("testWm_mc.root");
//TFile gg("../etmissb/Wmunu_noetcut.root");
//TFile gg("../etmissb/WenuMC_noetcut.root");
TFile gg("WmunuMC.root");

//TFile gg("pippo.root");
//TFile gg("../etmiss2/Wenu_noetcut_sumet.root");
//TFile gg("../etmissb/Wmunu_et25.root");
//TFile gg("Wmunu_noetcut_1_d3pd.root");
//TFile gg("testWm_mc.root");
TFile bz("ZmumuMC.root");

//TFile bb("result.root");
TFile bb("bbmu15_MC.root");
//TFile bw("../etmiss3/WtaunuMC.root");
//TFile bw("WtaunuinclMC_mu.root");
TFile bw("WtaunuMC_r1659.root");
TFile btt("ttbarMC.root");
TFile bWZ("WZMC.root");
TFile bWW("WWMC.root");
TTree* t=(TTree*)gg.Get("MiniT");

TFile fd("../etmissb/sumetratio.root");
TH1F *ww = (TH1F*)fd.Get("ratio");
float cutEt = 40.;
float cutMt = 0.;
float xcg[] = {cutEt,100.,100.,cutEt,cutEt};
float ycg[] = {120.,120.,cutMt,cutMt,120.};
TCutG etcut("etcut",5,xcg,ycg);
int metbin = 1;
//#define DATA



#ifndef __MAKECINT__
#ifdef DATA
TH1D *tofit =  ((TH2*)ff.Get("EtvsMt"))->ProjectionY("tofit",  metbin,100);
#else
TH1D *tofit =  ((TH2*)gg.Get("EtvsMt"))->ProjectionY("tofit",metbin,100);
#endif
TH1 *bkgz =   ((TH2*)bz.Get("EtvsMt"))->ProjectionY("bkgz",metbin,100);
TH1 *bkgb =   ((TH2*)bb.Get("EtvsMt"))->ProjectionY("bkgb",metbin,100);
//TH1 *bkgb =   ((TH2*)bbnoiso.Get("noniso"))->ProjectionY("bkgb",metbin,100);
TH1 *bkgw =  ((TH2*)bw.Get("EtvsMt"))->ProjectionY("bkgw", metbin,100);
TH1 *bkgtt =  ((TH2*)btt.Get("EtvsMt"))->ProjectionY("bkgtt", metbin,100);
TH1 *bkgWW =  ((TH2*)bWW.Get("EtvsMt"))->ProjectionY("bkgWW", metbin,100);
TH1 *bkgWZ =  ((TH2*)bWZ.Get("EtvsMt"))->ProjectionY("bkgWZ", metbin,100);
#else
TH1D *tofit,*bkgz,*bkgw,*bkgb,*bkgtt,*bkgWW,*bkgWZ;
#endif



//TH1F* tofit=(TH1F*)gg.Get("Mtref");
//TH2F* reso=(TH2F*)gg.Get("EtvsEtTrue");
TH1F* mcfit;
TH1F* mcfitd,zmm,wtn,bbm;
TH1F* mcfitsb;
TH2F* maxwell;

string pnames[] = { "scale","reso","bkgz","bkgb"};
double pmin[] = {0.1,0.01,0.0,0.0};
double pini[] = {1.00,0.52,0.05,0.1};
double pmax[] = {10.1,5. ,10.0,20.0};
double ss1[3000000];
double ss2[3000000];
double ss3[3000000];
int isme=0;


void fcn(int &npar, double *gin, double &f, double *par, int iflag);

void loasme(){
 for (int i = 0; i <3000000 ; i++) {
   ss1[i]=gRandom->Gaus(0,1);
   ss2[i]=gRandom->Gaus(0,1);
   ss3[i]=gRandom->Gaus(0,1);
 }
}
void minit(int option=1, const bool& mc=false){

  //The default minimizer is Minuit, you can also try Minuit2
  loasme();
  if (option%10 == 2) 
    TVirtualFitter::SetDefaultFitter("Minuit2");
  else
    TVirtualFitter::SetDefaultFitter("Minuit");

  TVirtualFitter * minuit = TVirtualFitter::Fitter(0,5);
  // set number of  parameters
  for (int i = 0; i < Npar; i++) {
    minuit->SetParameter(i, pnames[i].c_str(), pini[i], 1e-2*(pmax[i]-pmin[i]), pmin[i], pmax[i]);
  }
  minuit->SetFCN(fcn);
  //minuit->FixParameter(0);
  //         minuit->FixParameter(1);
#ifndef DATA
  minuit->FixParameter(2);
  minuit->FixParameter(3);
#endif
  //  minuit->FixParameter(4);
    
  double arglist[100];
  arglist[0] = 0;
  // set print level
  minuit->ExecuteCommand("SET PRINT",arglist,2);
  arglist[0] = 0.005;
  //  arglist[0] = 0.2;
  // set eps level
  minuit->ExecuteCommand("SET EPS",arglist,2);

  // minimize
  arglist[0] = 5000; // number of function calls
  arglist[1] = 0.05; // tolerance

  minuit->ExecuteCommand("SIMPLEX",arglist,2);
  minuit->ExecuteCommand("MIGRAD",arglist,2);
  minuit->ExecuteCommand("HESSE",arglist,2);

  //  minuit->ExecuteCommand("MINOS",arglist,2);
// //  arglist[0] = 3;
  minuit->ExecuteCommand("CALL FCN",arglist,1);

  //get result
  double minParams[10];
  double parErrors[10];
  for (int i = 0; i < Npar; ++i) {
    minParams[i] = minuit->GetParameter(i);
    parErrors[i] = minuit->GetParError(i);
  }

  double chi2, edm, errdef;
  int nvpar, nparx;
  minuit->GetStats(chi2,edm,errdef,nvpar,nparx);

  TMatrixD cov(2,2,TVirtualFitter::GetFitter()->GetCovarianceMatrix() ); 
  cov.Print();
  /*
  for (int i=0; i<2; i++){
    for (int j=0; j<2; j++){
      
      cout<<i<<" "<<j<<" "<< minuit->GetCovarianceMatrixElement(i,j)<<endl;
      //sqrt(minuit->GetCovarianceMatrixElement(j,j))/sqrt(minuit->GetCovarianceMatrixElement(i,i))<<endl;
      }						 
    }
  */

  int minbin = 1;
  int maxbin = tofit->GetNbinsX();
  double totmc = mcfitsb->Integral(minbin,maxbin);
 double totdt = tofit->Integral(minbin,maxbin);

  //  bkgz->Scale(totdt/totmc*minParams[2]);
  //  bkgw->Scale(totdt/totmc*minParams[4]);
  //  bkgb->Scale(totdt/totmc*minParams[3]);

}
void fcn(int &npar, double *gin, double &f, double *par, int iflag){
  //calculate chisquare

  f = 0.;
  isme=0;
  mcfit = (TH1F*)tofit->Clone();
  mcfit->Reset();
  mcfitsb = (TH1F*)tofit->Clone();
  mcfitsb->Reset();
  //  maxwell = (TH2F*)reso->Clone();

  TTree *fChain=t;
  if (fChain == 0) return;
  Long64_t nentries = fChain->GetEntriesFast();
  //  cout<<"dentro il loop"<<nentries<<endl;
  int a;
  //  cin >> a ;
  Float_t  EtMCx,EtMCy,Etrecx,Etrecy,sumEt,pxlep,pylep,etcone40;
  fChain->SetBranchAddress("EtMCx",&EtMCx);
  fChain->SetBranchAddress("EtMCy",&EtMCy);
  fChain->SetBranchAddress("Etrecx",&Etrecx);
  fChain->SetBranchAddress("Etrecy",&Etrecy);
  fChain->SetBranchAddress("sumEt",&sumEt);
  fChain->SetBranchAddress("pxlep",&pxlep);
  fChain->SetBranchAddress("pylep",&pylep);
  fChain->SetBranchAddress("etcone40",&etcone40);
  //  cout <<"par"<<" "<<par[0]<<endl;
  Long64_t nbytes = 0, nb = 0;
  //  int nentro = 135000;
  for (Long64_t jentry=0; jentry<nentries;jentry++) {
    
    //for (Long64_t jentry=1; jentry<nentro;jentry++) {
    Long64_t ientry = fChain->GetEntry(jentry);
    
    if (ientry < 0) break;
    nb = fChain->GetEntry(jentry);   nbytes += nb;
    
    //    cout<<jentry<<endl;

    double r = par[0];
    
    //    double r2 = par[2];
    double sme = fabs(par[1]);
       
    int ntest=5;
    
    double etxs,etys;   
    double etxsn,etysn;   


    TLorentzVector ETlep,ptot,ETMiss;
    TLorentzVector ptotrec,ETMissrec;
    ETlep.SetXYZM(pxlep,pylep,0,0);

    for(int kk=0; kk<ntest;kk++) {
      double scala=1.;
      //      etxs = r*EtMCx/fabs(EtMCx)+EtMCx+ss1[isme]*sme*pow(fabs(sumEt)*scala,0.5);
      //      etys = r*EtMCy/fabs(EtMCy)+EtMCy+ss2[isme]*sme*pow(fabs(sumEt)*scala,0.5);
      
      etxs = r*EtMCx +ss1[isme]*sme*pow(fabs(sumEt+ETlep.Et())*scala,0.5);
      etys = r*EtMCy+ ss2[isme]*sme*pow(fabs(sumEt+ETlep.Et())*scala,0.5);

      etxsn = EtMCx+ss1[isme]*sme*pow(fabs(sumEt)*scala,0.5);
      etysn = EtMCy+ss2[isme]*sme*pow(fabs(sumEt)*scala,0.5);
      
      /*      etxs = r*EtMCx+ss1[isme]*sme*pow(fabs(EtMCx),0.5);
      etys = r*EtMCy+ss2[isme]*sme*pow(fabs(EtMCy),0.5);

      etxsn = EtMCx+ss1[isme]*sme*pow(fabs(EtMCx),0.5);
      etysn = EtMCy+ss2[isme]*sme*pow(fabs(EtMCy),0.5);
      */
      ETMiss.SetXYZM(etxs,etys,0,0);
      ETMissrec.SetXYZM(Etrecx,Etrecy,0,0);
      ptot =  ETlep + ETMiss;
      ptotrec =  ETlep + ETMissrec;


      //      etxs = r*(EtMCx+ss1[isme]*sme*pow(fabs(EtMCx),0.5));
      //      etys = r*(EtMCy+ss2[isme]*sme*pow(fabs(EtMCy),0.5));
      //      cout <<EtMCx<<"->x "<<etxs <<endl;
      //      cout <<EtMCy<<"->y "<<etys <<endl;
      //      cout <<r<<"   "<<sme<<endl;
      //      cout <<ss1[isme]<<"   "<<ss2[isme]<<endl;
      //      cout <<"-----------------------"<<endl;

      isme += 1; 
      //      if(isme>=ntest*nentro)isme=0;         
      if(isme>=ntest*nentries)isme=0;         

      double etmiss = pow(etxs*etxs+etys*etys,0.5);
      double etmissn = pow(etxsn*etxsn+etysn*etysn,0.5);
      double etmissrec = pow(Etrecx*Etrecx+Etrecy*Etrecy,0.5);
      //      if(etmissrec<(metbin-1)*2.0)continue;
      if(etmiss<(metbin-1)*2.0)continue;

      //      if(etmissrec<xcg[0])continue;
      //    if (ptot.Mag() < ycg[2]) continue;
      //      cout<<ptot.Mag()<<" "<<ptotrec.Mag()<<endl;
      //      if(etmissrec<25)continue;

     
      double etmiss_true = pow(EtMCx*EtMCx+EtMCy*EtMCy,0.5);
      //      maxwell->Fill(etmiss_true,(etmissn-etmiss_true)/etmiss_true);



      int wwbin = int(sumEt/ww->GetBinWidth(1));
      
      if (wwbin > ww->GetNbinsX()) wwbin =  ww->GetNbinsX();
      double weighta = ww->GetBinContent(wwbin);
      double weighte = ww->GetBinError(wwbin);
      double  weight =weighta+ss3[isme]*weighte;
      //      cout<<"weight "<<weight<<endl;
      //      weight = 1.;
      //      weight =weighta;
#ifdef DATA
      mcfit->Fill(ptot.Mag(),1./ntest*weight*10);     
#else
      mcfit->Fill(ptot.Mag(),1./ntest);     
#endif
      //mcfit->Fill(ptot.Mag(),1./ntest);     
      //mcfit->Fill(etmiss,1./ntest);     
    }
  }
  //  int minbin = 11;
  int minbin = 11;
  int maxbin = tofit->GetNbinsX();
  

  cout<<"maxbin"<<maxbin<<endl;   
  //  cin >> a ;
  for(int j=minbin; j<=maxbin; j++){
    double bkg=0.;
    //    if(j>=20&&j<=100){ bkg= (-(j+20)+100.)*mcfit->Integral(minbin,maxbin)/tofit->Integral(minbin,maxbin); } 
    bkg=0.;

    double normttbar = 1.4413*1e-1*0.558/99928;
    double normWW = 2.9592*1e-2*0.38863/247868;
    double normWZ = 1.123*1e-2*0.30847/59971;
    double normWtau = 7.8/199975;
    //    double normZmm = 0.86/99989.;
    double normZmm = 0.86/100985.;
    
    //    double normWoverZ = 100985./199975*7.8/.86;
    double normWoverZ = 99989./199975*7.8/.86;
    //    double normWinclZ = 99989./1996438.*8.9/.86;
    //    double bkg1 = (bkgz->GetBinContent(j)+bkgw->GetBinContent(j)*normWoverZ)*par[2]/(bkgz->Integral(minbin,maxbin)+bkgw->Integral(minbin,maxbin)*normWoverZ);

    double bkg1 = par[2]*(bkgz->GetBinContent(j)*normZmm + bkgw->GetBinContent(j)*normWtau + bkgtt->GetBinContent(j)*normttbar +  bkgWZ->GetBinContent(j)*normWZ + bkgWW->GetBinContent(j)*normWW)/(bkgz->Integral(minbin,maxbin)*normZmm + bkgw->Integral(minbin,maxbin)*normWtau + bkgtt->Integral(minbin,maxbin)*normttbar +  bkgWZ->Integral(minbin,maxbin)*normWZ + bkgWW->Integral(minbin,maxbin)*normWW); 
      
    double bkg2 = bkgb->GetBinContent(j)*par[3]/bkgb->Integral(minbin,maxbin);
    

    //    cout<<bkg2 <<" "<<par[3]<<" "<<bkgb->GetBinContent(j)<<" "<<bkgb->Integral(minbin,maxbin)<<endl;
#ifdef DATA
    mcfitsb->SetBinContent(j, (mcfit->GetBinContent(j)/mcfit->Integral(minbin,maxbin)+bkg1+bkg2)*mcfit->Integral(minbin,maxbin));
#else
    mcfitsb->SetBinContent(j, mcfit->GetBinContent(j));
#endif
  } 

  for(int j=minbin; j<=maxbin; j++){
    double zbin  = tofit->GetBinContent(j)/tofit->Integral(minbin,maxbin);
    double ezbin = tofit->GetBinError(j)/tofit->Integral(minbin,maxbin);  
    //    if(j>=40&&j<=100){
    //      fbin+= (-j+100.)/mcfit->Integral(minbin,maxbin);}      
    double fbin  = mcfitsb->GetBinContent(j)/mcfitsb->Integral(minbin,maxbin);
    double efbin = mcfitsb->GetBinError(j)/mcfitsb->Integral(minbin,maxbin);
    
    //    cout<<fbin<<" bin "<<zbin<<endl;
    //    cout<<efbin<<" ebin "<<ezbin<<endl;


    if(tofit->GetBinContent(j)>20){
      double term = pow(fbin - zbin,2)/(ezbin*ezbin + efbin*efbin);	
      //      cout<<" term "<<term<<endl;
      f += term;
    }
  }
  double totmc = mcfitsb->Integral(minbin,maxbin);
  double totdt = tofit->Integral(minbin,maxbin);
  mcfitd = (TH1F*)mcfitsb->Clone();
  mcfitd->Scale(totdt/totmc);
  
 
  cout<<" chi2 "<<f<<" r "<<par[0]<<" reso "<<par[1]<<endl<<flush;
  //  cin >> a ;
}