#include "TDSubAna.hh"

int TDSubAna::HistFill(const char * histname, double fillvalue, double weight)
{
  //if(diltype==17)  //used for test
    /*if(_subAnaHists->GetAnaHist(histname)!=_subAnaHists->GetAnaHistFast(histname))
    {std::cout<<"Fast not working"<<std::endl;}
    else{std::cout<<"Fast working:"<<_subAnaHists->GetAnaHist(histname)->GetHistName()<<" fast: "<<_subAnaHists->GetAnaHistFast(histname)->GetHistName()<<std::endl;}*/
    /*TDAnaHist* histresult=_subAnaHists->GetAnaHistFast(histname);//speedup
    if(histresult!=NULL)
    {
      histresult->GetHist()->Fill(fillvalue, weight);
      return 0;
    }*/

  
    TDAnaHist* histresult=_subAnaHists->GetAnaHistFast(histname);
    if(histresult!=NULL) 
    {
      histresult->GetHist()->Fill(fillvalue, weight); 
      return 0;
    }
    std::cout<<"Histogram was not found to fill"<<std::endl;
    return 1;
}

int TDSubAna::HistFill(const char * histname, double fillvalue, double fillvalue2, double weight)
{
  //if(diltype==17)  //used for test
    if(_subAnaHists->GetAnaHist(histname)->GetHist()!=NULL)
    {
      ((TH2F*)(_subAnaHists->GetAnaHist(histname))->GetHist())->Fill(fillvalue, fillvalue2, weight);
      return 0;
    }
    std::cout<<"Histogram was not found to fill"<<std::endl;
    return 1;
}


bool TDSubAna::SameSign()
{    
      int charge=0;
      for(int lep=0; lep<nleptons; ++lep)
      {
        charge+=lepCharge[lep];
      }
      if(charge==0)return false; 
   return true;
}

double TDSubAna::GetDilWeight()
{
  LeptonType::lepton_type atype=LeptonType::UNDEFINED;
  LeptonType::lepton_type btype=LeptonType::UNDEFINED;
  double a_Et=0.0;
  double b_Et=0.0;
  bool firstlepfound=false;
  bool secondlepfound=false;
  double a_Eta=0.0;
  double b_Eta=0.0;
  for(int leps=0; leps < nleptons; ++leps)
  {
    if(secondlepfound)continue;
    //if(lepType[leps]==3 || lepType[leps]==0)
    
      if(!firstlepfound)
      {
        firstlepfound=true;
        atype=_leptonType.Int2Type(lepType[leps]);
	a_Et=lepEt[leps];
	a_Eta=lepEtaDet[leps];
      }else{
	secondlepfound=true;
	btype=_leptonType.Int2Type(lepType[leps]);
	b_Et=lepEt[leps];
	b_Eta=lepEtaDet[leps];
      } 
  }

  /*  //treat SCMIO as CMIO for now
      if(atype==LeptonType::SCMIO)atype=LeptonType::CMIO;
      if(btype==LeptonType::SCMIO)btype=LeptonType::CMIO;*/

  DilType::dilepton_type typedilType=_dilType.Int2Type(diltype);
  /*  return _dilWeight->DileptonWeight(_isdata, run,
      typedilType, atype, btype,
      a_Et, b_Et, 0.1, 0.1,
      metmag, a_Eta, b_Eta);*/
  /*  if(GetGoodGenZ()!=0 && GetGoodGenZ(diltype) !=0)            // information about events generated in the respective good runlists is available
      return (_dilWeight->DileptonWeight(_isdata, run,
      typedilType, atype, btype,
      a_Et, b_Et, 0.1, 0.1,
      metmag, a_Eta, b_Eta, true))*GetGoodGenZ()/GetGoodGenZ(diltype);*/
  return _dilWeight->DileptonWeight(_isdata, run,
      typedilType, atype, btype,
      a_Et, b_Et, 0.1, 0.1,
      metmag, a_Eta, b_Eta, true);

}

double TDSubAna::DeanDetEta(double const eta, double const z0)
{
  double detEta;
  double zprime = z0 + 184.15 / (tan(2 * atan(exp(-eta))));
  if (zprime >= 0) {
    detEta = -log(
	tan(
	  0.5 * atan(
	    184.15/(z0+184.15/tan(2*atan(exp(-eta))))
	    )
	  )
	);
  } else {
    detEta = -log(
	tan(
	  0.5 * (3.142 + atan(
	      184.15/(z0+184.15/tan(2*atan(exp(-eta))))
	      ))
	  )
	);
  }
  return detEta;
}

void TDSubAna::InitTMVA()
{
  // ---------------------------------------------------------------
  // choose MVA methods to be trained + tested
  Use_Cuts            = 0;
  Use_Likelihood      = 0;//1;
  Use_LikelihoodD     = 0;//1;
  Use_PDERS           = 0;//1;
  Use_HMatrix         = 0;//1;
  Use_Fisher          = 0;//1;
  Use_CFMlpANN        = 0;//1;
  Use_TMlpANN         = 1;//0;
  Use_BDT             = 0;//1;
  //Use_BDT_GiniIndex   = 0; // default BDT method
  //Use_BDT_CrossEntro  = 0;
  //Use_BDT_SdivStSpB   = 0;
  //Use_BDT_MisClass    = 0;
  //Use_BDT_Bagging_Gini= 0;
  // ---------------------------------------------------------------
  EvaluateVariables  = 0; // perform evaluation for each input variable
  // ---------------------------------------------------------------
  inputVars = new vector<std::string>;
  inputVars->push_back("LepAEn");
  inputVars->push_back("LepBEn");
  inputVars->push_back("dilmass");
  inputVars->push_back("metmag");
  //inputVars->push_back("ChargeProduct");
  inputVars->push_back("addEt");
  inputVars->push_back("dPhiMetLJ");
  inputVars->push_back("ntightjets");
  inputVars->push_back("jet1_Et");
  inputVars->push_back("jet2_Et");
  inputVars->push_back("dildphi");
  inputVars->push_back("lepR");
  inputVars->push_back("metsig");
  //     inputVars->push_back("diltype");

  reader_ee = new TMVA::Reader( *inputVars );
  reader2_ee = new TMVA::Reader( *inputVars );

  reader_emu = new TMVA::Reader( *inputVars );
  reader2_emu = new TMVA::Reader( *inputVars );

  reader_mumu = new TMVA::Reader( *inputVars );
  reader2_mumu = new TMVA::Reader( *inputVars );


  string dir    = "./weights/";
  char prefix[200];

  sprintf(prefix,"TMVTEST1_%d_0", 160);//only one training for now
  if (Use_TMlpANN)    reader_ee->BookMVA( "TMVA::Types::TMlpANN",      dir + prefix + "_TMlpANNBFGS.weights.txt");//"_TMlpANN.weights"  ); 
  sprintf(prefix,"TMVTEST2_%d_0", (int)_higgsMass);
  if (Use_TMlpANN)    reader2_ee->BookMVA( "TMVA::Types::TMlpANN",      dir + prefix + "_TMlpANNBFGS.weights.txt");//"_TMlpANN.weights"  );
  sprintf(prefix,"TMVTEST1_%d_1", 160);//only one training for now
  if (Use_TMlpANN)    reader_emu->BookMVA( "TMVA::Types::TMlpANN",      dir + prefix + "_TMlpANNBFGS.weights.txt");//"_TMlpANN.weights"  );        
  sprintf(prefix,"TMVTEST2_%d_1", (int)_higgsMass);        
  if (Use_TMlpANN)    reader2_emu->BookMVA( "TMVA::Types::TMlpANN",      dir + prefix + "_TMlpANNBFGS.weights.txt");//"_TMlpANN.weights"  );
  sprintf(prefix,"TMVTEST1_%d_2", 160);//only one training for now
  if (Use_TMlpANN)    reader_mumu->BookMVA( "TMVA::Types::TMlpANN",      dir + prefix + "_TMlpANNBFGS.weights.txt");//"_TMlpANN.weights"  );        
  sprintf(prefix,"TMVTEST2_%d_2", (int)_higgsMass);        
  if (Use_TMlpANN)    reader2_mumu->BookMVA( "TMVA::Types::TMlpANN",      dir + prefix + "_TMlpANNBFGS.weights.txt");//"_TMlpANN.weights"  );


  initTMVA=true;
}

bool TDSubAna::evaluateTMVA(double weight, TDSubAnaCutInterface * cutInterface, const char * TMLP2Hist, bool fillTMLPhist)
{
  //if (fillTMLPhist)
  //_hprof->Fill(*(cutInterface->dilmass()),*(cutInterface->LepBEn()),weight);

  //std::cout<<"in evaluate: dildphi: "<<dildphi<<std::endl;

  bool debugrun=true;
  vector<double> varValues;
  varValues.push_back(*(cutInterface->LepAEn()));
  varValues.push_back(*(cutInterface->LepBEn()));
  varValues.push_back(*(cutInterface->dilmass()));
  varValues.push_back(*(cutInterface->metmag()));
  varValues.push_back(*(cutInterface->addEt()));
  varValues.push_back(*(cutInterface->dPhiMetLJ()));
  varValues.push_back(*(cutInterface->ntightjets()));
  varValues.push_back(*(cutInterface->jet1_Et()));
  varValues.push_back(*(cutInterface->jet2_Et()));
  varValues.push_back(dildphi);
  varValues.push_back(_lepR);
  varValues.push_back(metsig);

  /*std::cout<<"NN values:"<<std::endl;
    if(debugrun)for(vector<double>::iterator it=varValues.begin(); it!=varValues.end(); it++)
    {
    std::cout<<*it<<std::endl;
    }*/

  double MVA_Likelihood; 
  double MVA_LikelihoodD;
  double MVA_PDERS;      
  double MVA_HMatrix;       
  double MVA_Fisher;       
  double MVA_CFMlpANN;      
  double MVA_TMlpANN;        
  double MVA_BDT;          

  if(_NNdiltype==0)if (Use_TMlpANN    )  MVA_TMlpANN     = reader_ee->EvaluateMVA( varValues, "TMVA::Types::TMlpANN"    );
  //  std::cout<<"here"<<std::endl;
  if(_NNdiltype==1)if (Use_TMlpANN    )  MVA_TMlpANN     = reader_emu->EvaluateMVA( varValues, "TMVA::Types::TMlpANN"    );
  //  std::cout<<"here"<<std::endl;  
  if(_NNdiltype==2)if (Use_TMlpANN    )  MVA_TMlpANN     = reader_mumu->EvaluateMVA( varValues, "TMVA::Types::TMlpANN"    );
  //  std::cout<<"here"<<std::endl;


  double MVA_Likelihood2;
  double MVA_LikelihoodD2;
  double MVA_PDERS2;
  double MVA_HMatrix2;
  double MVA_Fisher2;
  double MVA_CFMlpANN2;
  double MVA_TMlpANN2;
  double MVA_BDT2;

  if(_NNdiltype==0)if (Use_TMlpANN    )  MVA_TMlpANN2     = reader2_ee->EvaluateMVA( varValues, "TMVA::Types::TMlpANN"    );
  // std::cout<<"here"<<std::endl;

  if(_NNdiltype==1)if (Use_TMlpANN    )  MVA_TMlpANN2     = reader2_emu->EvaluateMVA( varValues, "TMVA::Types::TMlpANN"    );
  // std::cout<<"here"<<std::endl;

  if(_NNdiltype==2)if (Use_TMlpANN    )  MVA_TMlpANN2     = reader2_mumu->EvaluateMVA( varValues, "TMVA::Types::TMlpANN"    );
  // std::cout<<"here"<<std::endl;



  if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA_TMlpANN"    , MVA_TMlpANN    , weight);
  if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA_TMlpANN_fine",MVA_TMlpANN    , weight);
  if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA2_TMlpANN"    , MVA_TMlpANN2    , weight);

  if(_NNdiltype==0)if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA_TMlpANN_ee"    , MVA_TMlpANN    , weight);
  if(_NNdiltype==0)if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA_TMlpANN_fine_ee",MVA_TMlpANN    , weight);
  if(_NNdiltype==0)if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA2_TMlpANN_ee"    , MVA_TMlpANN2    , weight);
  if(_NNdiltype==1)if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA_TMlpANN_emu"    , MVA_TMlpANN    , weight);
  if(_NNdiltype==1)if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA_TMlpANN_fine_emu",MVA_TMlpANN    , weight);
  if(_NNdiltype==1)if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA2_TMlpANN_emu"    , MVA_TMlpANN2    , weight);
  if(_NNdiltype==2)if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA_TMlpANN_mumu"    , MVA_TMlpANN    , weight);
  if(_NNdiltype==2)if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA_TMlpANN_fine_mumu",MVA_TMlpANN    , weight);
  if(_NNdiltype==2)if (Use_TMlpANN    && fillTMLPhist) HistFill("MVA2_TMlpANN_mumu"    , MVA_TMlpANN2    , weight);

  //if(_NNdiltype==0)
  //{
    if(Use_TMlpANN && fillTMLPhist && MVA_TMlpANN> GetNNCutVal(1, (int)_NNdiltype))++_passNNDYcut[1];
    if(Use_TMlpANN && fillTMLPhist && MVA_TMlpANN> GetNNCutVal(2, (int)_NNdiltype))++_passNNDYcut[2];
    if(Use_TMlpANN && fillTMLPhist && MVA_TMlpANN> GetNNCutVal(3, (int)_NNdiltype))++_passNNDYcut[3];
  //}

  if(Use_TMlpANN && fillTMLPhist && MVA_TMlpANN<0.8 && MVA_TMlpANN>0.6 && _drawNNhists4)
  {
    _histsNNDYCutMid.Fill(weight);
//#include "AnaHistFillNNDYCutMid.cc"
  }
  //if(debugrun)std::cout<<"DY, WW NN score:"<<MVA_TMlpANN<<" "<<MVA_TMlpANN2<<" "<<std::endl; 
  if(
      (Use_TMlpANN && MVA_TMlpANN >GetNNCutVal(2, _NNdiltype)  && _7708Cutrejection==0) ||           //GetNNCutVal(2, _NNdiltype) or GetNNsdyrbCutVal(_NNdiltype) 
      (_7708Cutrejection>0 && _extraDYpass ) ||
      (Use_TMlpANN && _7708Cutrejection<0 && _extraDYpass && MVA_TMlpANN> GetNNCutVal(2, _NNdiltype))   //GetNNCutVal(2, _NNdiltype) or GetNNsdyrbCutVal(_NNdiltype)
    )
  {
    if (Use_TMlpANN    )                 HistFill(TMLP2Hist    , MVA_TMlpANN2    , weight);  // hist filling of the JES ans LES histograms   
    if (Use_TMlpANN    && fillTMLPhist ) HistFill("MVA_TMlpANN_cut1"    , MVA_TMlpANN    , weight);
    if (Use_TMlpANN    && fillTMLPhist )
    {
      //if(debugrun)std::cout<<"in anadump: "<<MVA_TMlpANN2<<" "<<MVA_TMlpANN<<" "<<weight<<" "<<GetScaleforNormalisation()<<std::endl;
      //if(_initFakes)_anafiledump<<run<<" "<<evt<<" "<<_NNdiltype<<" "<<_leptonType.Int2String(lepType[0])<<" "<<_leptonType.Int2String(lepType[1])<<" "<<MVA_TMlpANN2<<" "<<weight*GetScaleforNormalisation()<<" "<<weight*(1+GetRelativeFakeRateError(run, lepType[1], lepEt[1]))<<" "<<weight*(1-GetRelativeFakeRateError(run, lepType[1], lepEt[1]))<<std::endl;  //fakes have the fakeerror weights with them
      //else _anafiledump<<run<<" "<<evt<<" "<<_NNdiltype<<" "<<_leptonType.Int2String(lepType[0])<<" "<<_leptonType.Int2String(lepType[1])<<" "<<MVA_TMlpANN2<<" "<<weight*GetScaleforNormalisation()<<std::endl;
      float sexoweight=1.0;
      double shichehspecialmet=metmag*sin(dPhiMetLJ);
      if(_sexocombined && dilmass<25)for(int lep=0; lep<nleptons; ++lep)
      {
	sexoweight/=_sexoweights[lepType[lep]];
      }

      if(_initFakes)_anafiledump<<run<<" "
	                          <<evt<<" "
	                            <<_NNdiltype<<" "
	                            <<_leptonType.Int2String(lepType[0])<<" "
	                            <<_leptonType.Int2String(lepType[1])<<" "
	                            <<MVA_TMlpANN2<<" "
	                            <<weight<<" "
	                            <<1.0<<" "
	                            <<GetFakeRate(run, lepType[1], lepEt[1])<<" "
	                            <<(1+GetRelativeFakeRateError(run, lepType[1], lepEt[1]))<<" "
	                            <<(1-GetRelativeFakeRateError(run, lepType[1], lepEt[1]))<<" "
	                            <<GetXS()<<" "
	                            <<GetInputEvents(diltype)<<" "
	                            <<*(cutInterface->ntightjets())<<" "
	                            <<_alldilLums[diltype]<<" "
	                            <<rawweight<<" "
	                            <<sexoweight<<" "
	                            <<metmag<<" "
	                            <<metsig<<" "
	                            <<shichehspecialmet<<" "
	                            <<dilmass<<" "
	                            <<dildphi<<" "
	                            <<lepEt[0]<<" "
	                            <<lepEt[1]<<" "
	                            <<lepEtaDet[0]<<" "
	                            <<lepEtaDet[1]<<" "
	                            <<std::endl;  //fakes have the fakeerror weights with them
      else if(isData()) _anafiledump<<run<<" "
	                            <<evt<<" " 
	                            <<_NNdiltype<<" "
	                            <<_leptonType.Int2String(lepType[0])<<" "
	                            <<_leptonType.Int2String(lepType[1])<<" "
	                            <<MVA_TMlpANN2<<" "
	                            <<weight<<" " 
	                            <<1.0<<" "
	                            <<1.0<<" "
	                            <<1.0<<" "
	                            <<1.0<<" "
	                            <<1.0<<" "
	                            <<1.0<<" "
	                            <<*(cutInterface->ntightjets())<<" "
	                            <<_alldilLums[diltype]<<" "
	                            <<rawweight<<" "
	                            <<sexoweight<<" "
	                            <<metmag<<" "
	                            <<metsig<<" "
	                            <<shichehspecialmet<<" "
	                            <<dilmass<<" "
	                            <<dildphi<<" "
	                            <<lepEt[0]<<" "
	                            <<lepEt[1]<<" "
	                            <<lepEtaDet[0]<<" "
	                            <<lepEtaDet[1]<<" "
	                            <<std::endl;
      else              _anafiledump<<run<<" "
	                          <<evt<<" " 
	                            <<_NNdiltype<<" "
	                            <<_leptonType.Int2String(lepType[0])<<" "
	                            <<_leptonType.Int2String(lepType[1])<<" "
	                            <<MVA_TMlpANN2<<" "
	                            <<weight<<" " 
	                            <<GetScaleforNormalisation()<<" "
	                            <<1.0<<" "
	                            <<1.0<<" "
	                            <<1.0<<" "
	                            <<GetXS()<<" "
	                            <<GetInputEvents(diltype)<<" "
	                            <<*(cutInterface->ntightjets())<<" "
	                            <<_alldilLums[diltype]<<" "
	                            <<rawweight<<" "
	                            <<sexoweight<<" "
	                            <<metmag<<" "
	                            <<metsig<<" "
	                            <<shichehspecialmet<<" "
	                            <<dilmass<<" "
	                            <<dildphi<<" "
	                            <<lepEt[0]<<" "
	                            <<lepEt[1]<<" "
	                            <<lepEtaDet[0]<<" "
	                            <<lepEtaDet[1]<<" "
                                    <<std::endl;
    }
    if(fillTMLPhist && _drawNNhists)
    {
      if(dilmass<120)
      {
	HistFill("HighMllControlNN_MVA_TMlpANN"    , MVA_TMlpANN    , weight);
	HistFill("HighMllControlNN_MVA2_TMlpANN"    , MVA_TMlpANN2    , weight);
	HistFill("HighMllControlNN_MVA2_TMlpANN_cut1"    , MVA_TMlpANN2    , weight);
      }
      if(dilmass>66 && dilmass<116)
      {
	HistFill("ZPeakNN_MVA_TMlpANN"    , MVA_TMlpANN    , weight);
	HistFill("ZPeakNN_MVA2_TMlpANN"    , MVA_TMlpANN2    , weight);
	HistFill("ZPeakNN_MVA2_TMlpANN_cut1"    , MVA_TMlpANN2    , weight);
      }
    }
    if(fillTMLPhist && _drawNNhists3)
    {
      if(MVA_TMlpANN2<0.1)
      {
	_histsWWNNlowCut.Fill(weight);
//#include "AnaHistFillWWNNlowCut.cc"
      }else if(MVA_TMlpANN2<0.5){
	_histsWWNNmidCut.Fill(weight);
//#include "AnaHistFillWWNNmidCut.cc"
      }else{
	 _histsWWNNhighCut.Fill(weight);
//#include "AnaHistFillWWNNhighCut.cc"
      }
    }
    if(fillTMLPhist && _drawNNhists2)
    {
      if(*(cutInterface->LepAEt())>60 && *(cutInterface->LepBEt())>50)
      {
	HistFill("LepEtControlNN_MVA_TMlpANN"    , MVA_TMlpANN    , weight);
        HistFill("LepEtControlNN_MVA2_TMlpANN"    , MVA_TMlpANN2    , weight);
        HistFill("LepEtControlNN_MVA2_TMlpANN_cut1"    , MVA_TMlpANN2    , weight);
      }
    }
    
    if(_NNdiltype==0)
    {
      if (Use_TMlpANN    && fillTMLPhist ) HistFill("MVA_TMlpANN_cut1_ee"    , MVA_TMlpANN    , weight);
      if (Use_TMlpANN    && fillTMLPhist ) HistFill("MVA2_TMlpANN_cut1_ee"    , MVA_TMlpANN2    , weight);
      if (Use_TMlpANN    && fillTMLPhist )
      {
	for(int pdf=0; pdf!=40; ++pdf)HistFill("MVA2_TMlpANN_cut1_PDF_ee", MVA_TMlpANN2, (float)pdf+0.5, weight*PDFWeight[5]/PDFWeight[pdf+6]);
	if(_initFakes){
	  HistFill("MVA2_TMlpANN_cut1_ISRup_ee"   , MVA_TMlpANN2, weight);
	  HistFill("MVA2_TMlpANN_cut1_ISRdn_ee"   , MVA_TMlpANN2, weight);
	  HistFill("MVA2_TMlpANN_cut1_Fakeup_ee"  , MVA_TMlpANN2, weight*(1+GetRelativeFakeRateError(run, lepType[1], lepEt[1])));
	  HistFill("MVA2_TMlpANN_cut1_Fakedn_ee"  , MVA_TMlpANN2, weight*(1-GetRelativeFakeRateError(run, lepType[1], lepEt[1])));
	}else{
	  HistFill("MVA2_TMlpANN_cut1_ISRup_ee"   , MVA_TMlpANN2, weight*(1+GetISRError(true, subprocesspT)));
	  HistFill("MVA2_TMlpANN_cut1_ISRdn_ee"   , MVA_TMlpANN2, weight*(1+GetISRError(false, subprocesspT)));
	  /*if(_datasource=="wgamma" && _wgammaerr)//now included as a np
	  {
	    HistFill("MVA2_TMlpANN_cut1_Fakeup_ee", MVA_TMlpANN2, weight*1.2);
	    HistFill("MVA2_TMlpANN_cut1_Fakedn_ee", MVA_TMlpANN2, weight*0.8);
	  }else{*/
	    HistFill("MVA2_TMlpANN_cut1_Fakeup_ee", MVA_TMlpANN2, weight);
	    HistFill("MVA2_TMlpANN_cut1_Fakedn_ee", MVA_TMlpANN2, weight);
	 // }
	}
      }
    }else if(_NNdiltype==1){
      if (Use_TMlpANN    && fillTMLPhist ) HistFill("MVA_TMlpANN_cut1_emu"    , MVA_TMlpANN    , weight);
      if (Use_TMlpANN    && fillTMLPhist ) HistFill("MVA2_TMlpANN_cut1_emu"    , MVA_TMlpANN2    , weight);
      if (Use_TMlpANN    && fillTMLPhist )
      {
	for(int pdf=0; pdf!=40; ++pdf)HistFill("MVA2_TMlpANN_cut1_PDF_emu", MVA_TMlpANN2, (float)pdf+0.5, weight*PDFWeight[5]/PDFWeight[pdf+6]);
	if(_initFakes){
	  HistFill("MVA2_TMlpANN_cut1_ISRup_emu"   , MVA_TMlpANN2, weight);
	  HistFill("MVA2_TMlpANN_cut1_ISRdn_emu"   , MVA_TMlpANN2, weight);
	  HistFill("MVA2_TMlpANN_cut1_Fakeup_emu"  , MVA_TMlpANN2, weight*(1+GetRelativeFakeRateError(run, lepType[1], lepEt[1])));
	  HistFill("MVA2_TMlpANN_cut1_Fakedn_emu"  , MVA_TMlpANN2, weight*(1-GetRelativeFakeRateError(run, lepType[1], lepEt[1])));
	}else{
	  HistFill("MVA2_TMlpANN_cut1_ISRup_emu"   , MVA_TMlpANN2, weight*(1+GetISRError(true, subprocesspT)));
	  HistFill("MVA2_TMlpANN_cut1_ISRdn_emu"   , MVA_TMlpANN2, weight*(1+GetISRError(false, subprocesspT)));
	  /*if(_datasource=="wgamma" && _wgammaerr)//now included as a np
	  {
	    HistFill("MVA2_TMlpANN_cut1_Fakeup_emu", MVA_TMlpANN2, weight*1.2);
	    HistFill("MVA2_TMlpANN_cut1_Fakedn_emu", MVA_TMlpANN2, weight*0.8);
	  }else{*/
	    HistFill("MVA2_TMlpANN_cut1_Fakeup_emu", MVA_TMlpANN2, weight);
	    HistFill("MVA2_TMlpANN_cut1_Fakedn_emu", MVA_TMlpANN2, weight);
	  //}
	}
      }
    }else if(_NNdiltype==2){
      if (Use_TMlpANN    && fillTMLPhist ) HistFill("MVA_TMlpANN_cut1_mumu"    , MVA_TMlpANN    , weight);
      if (Use_TMlpANN    && fillTMLPhist ) HistFill("MVA2_TMlpANN_cut1_mumu"    , MVA_TMlpANN2    , weight);
      if (Use_TMlpANN    && fillTMLPhist )
      {
	for(int pdf=0; pdf!=40; ++pdf)HistFill("MVA2_TMlpANN_cut1_PDF_mumu", MVA_TMlpANN2, (float)pdf+0.5, weight*PDFWeight[5]/PDFWeight[pdf+6]);
	if(_initFakes){
	  HistFill("MVA2_TMlpANN_cut1_ISRup_mumu"   , MVA_TMlpANN2, weight);
	  HistFill("MVA2_TMlpANN_cut1_ISRdn_mumu"   , MVA_TMlpANN2, weight);
	  HistFill("MVA2_TMlpANN_cut1_Fakeup_mumu"  , MVA_TMlpANN2, weight*(1+GetRelativeFakeRateError(run, lepType[1], lepEt[1])));
	  HistFill("MVA2_TMlpANN_cut1_Fakedn_mumu"  , MVA_TMlpANN2, weight*(1-GetRelativeFakeRateError(run, lepType[1], lepEt[1])));
	}else{
	  HistFill("MVA2_TMlpANN_cut1_ISRup_mumu"   , MVA_TMlpANN2, weight*(1+GetISRError(true, subprocesspT)));
	  HistFill("MVA2_TMlpANN_cut1_ISRdn_mumu"   , MVA_TMlpANN2, weight*(1+GetISRError(false, subprocesspT)));
	  /*if(_datasource=="wgamma" && _wgammaerr)//now included as a np
	  {
	    HistFill("MVA2_TMlpANN_cut1_Fakeup_mumu", MVA_TMlpANN2, weight*1.2);
	    HistFill("MVA2_TMlpANN_cut1_Fakedn_mumu", MVA_TMlpANN2, weight*0.8);
	  }else{*/
	    HistFill("MVA2_TMlpANN_cut1_Fakeup_mumu", MVA_TMlpANN2, weight);
	    HistFill("MVA2_TMlpANN_cut1_Fakedn_mumu", MVA_TMlpANN2, weight);
	  //}
	}
      }
    }
    return true;
  }
  return false;
}

void TDSubAna::InitTMLP()
{
  TTree * totaltree;
  char buffer[200];
  //_myNN=new (TMultiLayerPerceptron)("LepAEn,LepBEn,dilmass,metmag,ChargeProduct,addEt,dPhiMetLJ,ntightjets,jet1_Et,jet2_Et,dildphi,metsig,lepR,diltype:14:type","weight",totaltree,"Entry$%2","((Entry$+1)%2)");
    printf(buffer,"./results/TMLP_%d.wgt", (int)_higgsMass);
  //_myNN->LoadWeights(buffer);
  _initTMLP=true;
}
double TDSubAna::EvaluateTMLP(TDSubAnaCutInterface * cutInterface)
{
  Double_t params[13];
  params[0]=(*(cutInterface->LepAEn()));
  params[1]=(*(cutInterface->LepBEn()));
  params[2]=(*(cutInterface->dilmass()));
  params[3]=(*(cutInterface->metmag()));
  params[4]=(*(cutInterface->addEt()));
  params[5]=(*(cutInterface->dPhiMetLJ()));
  params[6]=(*(cutInterface->ntightjets()));
  params[7]=(*(cutInterface->jet1_Et()));
  params[8]=(*(cutInterface->jet2_Et()));
  params[9]=(dildphi);
  params[10]=(_lepR);
  params[11]=(metsig);
  params[12]=(_NNdiltype);
  return _myNN->Evaluate(0,params);
}

double TDSubAna::GetPDFError(bool high)//TD high=true for upper limit
{
  double _pdfset[46];//TD dummy array for compliler - to be relpaced with mtuple values containing the PDF weights.
  double sumsq=0;
  double max1=0;
  for(int set=6; set!=46; set+=2)
  {
    if(high) max1=max( (_pdfset[set]-_pdfset[5]), (_pdfset[set+1]-_pdfset[5]) );
    if(!high)max1=max( (_pdfset[5]-_pdfset[set]), (_pdfset[5]-_pdfset[set+1]) );
    if(max1>0)sumsq+=max1*max1;
  }
  return sqrt(sumsq); 
}

void TDSubAna::InitISRHist(const char * filename)
{
  TFile * ISRFile=new (TFile)(filename);
  _ISRHist = (TH1F*)ISRFile->Get("hr_hcen")->Clone(); _ISRHist->SetDirectory(0);
  _ISRHistUp = (TH1F*)ISRFile->Get("hr_h64up")->Clone(); _ISRHistUp->SetDirectory(0);
  _ISRHistDn = (TH1F*)ISRFile->Get("hr_h64down")->Clone(); _ISRHistDn->SetDirectory(0);
  CheckISRHist();
  ISRFile->Close();
}

void TDSubAna::CheckISRHist()
{
  std::cout<<"_ISRHist"<<_ISRHist->Integral()<<std::endl;
  std::cout<<"_ISRHistUp"<<_ISRHistUp->Integral()<<std::endl;
  std::cout<<"_ISRHistDn"<<_ISRHistDn->Integral()<<std::endl;
}

double TDSubAna::GetISRError(bool high, double pT)//TD high=true for upper limit
{
  if(_debug)std::cout<<"in get ISRError "<<high<<" "<<pT<<std::endl;
  double midval=_ISRHist->GetBinContent(_ISRHist->FindBin(pT));
  double errorval=0;
  
  if(high)errorval=_ISRHistUp->GetBinContent(_ISRHist->FindBin(pT));
  else    errorval=_ISRHistDn->GetBinContent(_ISRHist->FindBin(pT));
  
  //return fabs((errorval-midval)/midval);
  return (errorval-midval)/midval;
}

void TDSubAna::InitFakesHists(const char * fname)
{
  _initFakes=true;
  _fakevectord=new std::vector<TH1F*>();
  _fakevectorh=new std::vector<TH1F*>();
  _fakevectori=new std::vector<TH1F*>();
  
  TFile fRootFile(fname, "read");
  if (!fRootFile.IsOpen()) {
    std::cout << "TFakeRate: ERROR: Could not open file " << fname << std::endl;
    std::cout << "TFakeRate: All GetFakeRate() calls will now return -1" << std::endl;
  }
  std::cout<<"got here"<<std::endl;
//TD the lepton types: TCE,  PHX,  PEM,  CMUP,  CMX,  CMU,  CMP,  CMIO
//TD the lepton types: TCE,  PHX,  PEM,  CMUP,  CMX,  CMU,  CMP,  CMIO, BMU,  SCMIO, DSCMIO, DCMIO, UNDEFINED
  fRootFile.ls();
  if(  (TH1F*)(fRootFile.Get("FakeRates/CEM_d"))  ==NULL)std::cout<<"eh?"<<std::endl;
  _fakevectord->push_back(   (TH1F*)(fRootFile.Get("FakeRates/CEM_d"))->Clone()    );
  std::cout<<"got here"<<std::endl;

  _fakevectord->push_back((TH1F*)fRootFile.Get("FakeRates/PHX_d")->Clone());
  _fakevectord->push_back(NULL);//TD PEM - should never be called!!!
  _fakevectord->push_back((TH1F*)fRootFile.Get("FakeRates/CMUP_d")->Clone());
  _fakevectord->push_back((TH1F*)fRootFile.Get("FakeRates/CMX_d")->Clone());
  _fakevectord->push_back((TH1F*)fRootFile.Get("FakeRates/CMU_d")->Clone());
  _fakevectord->push_back((TH1F*)fRootFile.Get("FakeRates/CMP_d")->Clone());
  _fakevectord->push_back((TH1F*)fRootFile.Get("FakeRates/CMIO_d")->Clone());
  //std::cout<<"fake debug"<<std::endl;
  _fakevectord->push_back(NULL);//TD BMU - should never be called!!!
  //std::cout<<"fake debug"<<std::endl;
  _fakevectord->push_back(new (TH1F)(*(*_fakevectord)[7])); //SCMIO uses the same fake rate as CMIO.
  //std::cout<<"fake debug"<<std::endl;

  _fakevectorh->push_back((TH1F*)fRootFile.Get("FakeRates/CEM_h")->Clone());
  _fakevectorh->push_back((TH1F*)fRootFile.Get("FakeRates/PHX_h")->Clone());
  _fakevectorh->push_back(NULL);//TD PEM - shoulh never be calleh!!!
  _fakevectorh->push_back((TH1F*)fRootFile.Get("FakeRates/CMUP_h")->Clone());
  _fakevectorh->push_back((TH1F*)fRootFile.Get("FakeRates/CMX_h")->Clone());
  _fakevectorh->push_back((TH1F*)fRootFile.Get("FakeRates/CMU_h")->Clone());
  _fakevectorh->push_back((TH1F*)fRootFile.Get("FakeRates/CMP_h")->Clone());
  _fakevectorh->push_back((TH1F*)fRootFile.Get("FakeRates/CMIO_h")->Clone());
  _fakevectorh->push_back(NULL);//TD BMU - should never be called!!!
  _fakevectorh->push_back(new (TH1F)(*(*_fakevectorh)[7])); //SCMIO uses the same fake rate as CMIO.
    
 
  _fakevectori->push_back((TH1F*)fRootFile.Get("FakeRates/CEM_i")->Clone());
  _fakevectori->push_back((TH1F*)fRootFile.Get("FakeRates/PHX_i")->Clone());
  _fakevectori->push_back(NULL);//TD PEM - shouli never be callei!!!
  _fakevectori->push_back((TH1F*)fRootFile.Get("FakeRates/CMUP_i")->Clone());
  _fakevectori->push_back((TH1F*)fRootFile.Get("FakeRates/CMX_i")->Clone());
  _fakevectori->push_back((TH1F*)fRootFile.Get("FakeRates/CMU_i")->Clone());
  _fakevectori->push_back((TH1F*)fRootFile.Get("FakeRates/CMP_i")->Clone());
  _fakevectori->push_back((TH1F*)fRootFile.Get("FakeRates/CMIO_i")->Clone());
  _fakevectori->push_back(NULL);//TD BMU - should never be called!!!
  _fakevectori->push_back(new (TH1F)(*(*_fakevectori)[7])); //SCMIO uses the same fake rate as CMIO.
    
  std::cout<<"fake debug"<<std::endl;
  for(int i=0; i!=10; i++)
  {
    if(i!=2 && i!=8)//nopem no bmu
    { 
      if((*_fakevectord)[i]==NULL)std::cout<<"nullhisto"<<std::endl;
      std::cout<<"fake debug"<<i<<std::endl;
      (*_fakevectord)[i]->SetDirectory(0);
      (*_fakevectorh)[i]->SetDirectory(0);
      (*_fakevectori)[i]->SetDirectory(0);
    }
  }
  std::cout<<"fake debug"<<std::endl;  
  //if(_debug)
  for(int i=0; i!=10; i++)
  {
    if(i!=2 && i!=8)//nopem no bmu
    {
      std::cout<<(*_fakevectord)[i]<<std::endl;
      std::cout<<((*_fakevectord)[i])->Integral()<<std::endl;
      std::cout<<(*_fakevectorh)[i]<<std::endl;
      std::cout<<((*_fakevectorh)[i])->Integral()<<std::endl;
      std::cout<<(*_fakevectori)[i]<<std::endl;
      std::cout<<((*_fakevectori)[i])->Integral()<<std::endl;
    }
  }
  fRootFile.Close();
}

double TDSubAna::GetFakeRate(int run, int leptontype, double Et)
{
  DebugOut("TDSubAna::GetFakeRate");
  /*if(_debug)
    for(int i=0; i!=8; i++)
    {
      if(i!=2)//nopem
      {
	std::cout<<(*_fakevectord)[i]<<std::endl;
	std::cout<<((*_fakevectord)[i])->Integral()<<std::endl;
	std::cout<<(*_fakevectorh)[i]<<std::endl;
	std::cout<<((*_fakevectorh)[i])->Integral()<<std::endl;
	std::cout<<(*_fakevectori)[i]<<std::endl;
	std::cout<<((*_fakevectori)[i])->Integral()<<std::endl;
      }
    }*/
  if(_debug)std::cout<<" "<<run<<" "<<leptontype<<" "<<Et<<std::endl;
  if(run<186599){return ((*_fakevectord)[leptontype])->GetBinContent(((*_fakevectord)[leptontype])->FindBin(Et));}
  else if(run>203799){return ((*_fakevectori)[leptontype])->GetBinContent(((*_fakevectori)[leptontype])->FindBin(Et));}
  else{return ((*_fakevectorh)[leptontype])->GetBinContent(((*_fakevectorh)[leptontype])->FindBin(Et));}
  return 0.0;
}

double TDSubAna::GetFakeRateError(int run, int leptontype, double Et)
{
  if(_debug)std::cout<<"in get fake error "<<run<<" "<<leptontype<<" "<<Et<<std::endl;
  if(run<186599){return ((*_fakevectord)[leptontype])->GetBinError(((*_fakevectord)[leptontype])->FindBin(Et));}
  else if(run>203799){return ((*_fakevectorh)[leptontype])->GetBinError(((*_fakevectorh)[leptontype])->FindBin(Et));}
  else{return ((*_fakevectori)[leptontype])->GetBinError(((*_fakevectori)[leptontype])->FindBin(Et));}
  return 0.0;
}

double TDSubAna::GetRelativeFakeRateError(int run, int leptontype, double Et)
{
  if(_debug)std::cout<<"in get Relative fake error "<<run<<" "<<leptontype<<" "<<Et<<std::endl;
  if(run<186599)
  {
    return ((*_fakevectord)[leptontype])->GetBinError(((*_fakevectord)[leptontype])->FindBin(Et))/
    ((*_fakevectord)[leptontype])->GetBinContent(((*_fakevectord)[leptontype])->FindBin(Et));
  }
  else if(run>203799)
  {
    return ((*_fakevectorh)[leptontype])->GetBinError(((*_fakevectorh)[leptontype])->FindBin(Et))/
    ((*_fakevectorh)[leptontype])->GetBinContent(((*_fakevectorh)[leptontype])->FindBin(Et));
  }
  else{
    return ((*_fakevectori)[leptontype])->GetBinError(((*_fakevectori)[leptontype])->FindBin(Et))/
    ((*_fakevectori)[leptontype])->GetBinContent(((*_fakevectori)[leptontype])->FindBin(Et));
  }
  return 0.0;
}

double TDSubAna::GetTransMass(std::vector<TLorentzVector*> thevec)
{
  TLorentzVector vect(0.0, 0.0, 0.0, 0.0);
  for(int i=0; i!=(int)thevec.size(); i++)
  {
    vect+=*(thevec[i]);
  }
  return vect.M();
}

void TDSubAna::RedoAxisLables()
{
  _subAnaHists->GetAnaHist("NNVar_LepAEt")->SetAxisLables("Higher Lepton Et (GeV)","Events");
  _subAnaHists->GetAnaHist("NNVar_LepAEt_ee")->SetAxisLables("Higher Lepton Et (GeV)","Events");
  _subAnaHists->GetAnaHist("NNVar_LepAEt_emu")->SetAxisLables("Higher Lepton Et (GeV)","Events");
  _subAnaHists->GetAnaHist("NNVar_LepAEt_mumu")->SetAxisLables("Higher Lepton Et (GeV)","Events");
  _subAnaHists->GetAnaHist("NNVar_LepBEt")->SetAxisLables("Lower Lepton Et GeV","Events");
  _subAnaHists->GetAnaHist("NNVar_LepBEt_ee")->SetAxisLables("Lower Lepton Et GeV","Events");
  _subAnaHists->GetAnaHist("NNVar_LepBEt_emu")->SetAxisLables("Lower Lepton Et GeV","Events");
  _subAnaHists->GetAnaHist("NNVar_LepBEt_mumu")->SetAxisLables("Lower Lepton Et GeV","Events");
  _subAnaHists->GetAnaHist("NNVar_LepAEn")->SetAxisLables("Higher Lepton Energy (GeV)","Events");
  _subAnaHists->GetAnaHist("NNVar_LepBEn")->SetAxisLables("Lower Lepton Energy GeV","Events");
  _subAnaHists->GetAnaHist("NNVar_dilmass")->SetAxisLables("Dilepton Mass (GeV)","Events");
  _subAnaHists->GetAnaHist("NNVar_metmag")->SetAxisLables("Missing Et (GeV)","Events");
  _subAnaHists->GetAnaHist("NNVar_addEt")->SetAxisLables("Sum Et (GeV)","Events");
  _subAnaHists->GetAnaHist("NNVar_dPhiMetLJ")->SetAxisLables("Closest object to Missing Et (Rads)","Events");
  _subAnaHists->GetAnaHist("NNVar_ntightjets")->SetAxisLables("Jet Number","Events");
  _subAnaHists->GetAnaHist("NNVar_jet1_Et")->SetAxisLables("Highest Jet Energy (GeV)","Events");
  _subAnaHists->GetAnaHist("NNVar_jet2_Et")->SetAxisLables("Second Highest Jet Energy (GeV)","Events");
  _subAnaHists->GetAnaHist("NNVar_dildphi")->SetAxisLables("Lepton Angular Separation (Rads)","Events");
  _subAnaHists->GetAnaHist("NNVar_lepR")->SetAxisLables("Dilepton #DeltaR","Events");
  _subAnaHists->GetAnaHist("NNVar_metsig")->SetAxisLables("Missing Et Significance","Events");
  _subAnaHists->GetAnaHist("WWmetsig")->SetAxisLables("Missing Et significance","Events");
  _subAnaHists->GetAnaHist("WWdilmass")->SetAxisLables("Dilepton Mass (GeV)","Events");
  _subAnaHists->GetAnaHist("WWmetmag")->SetAxisLables("Missing Et (GeV)","Events");
  _subAnaHists->GetAnaHist("WWsumEt")->SetAxisLables("Total Et (GeV)","Events");
  _subAnaHists->GetAnaHist("WWlepPt")->SetAxisLables("Lepton pt (GeV)","# leptons");
  _subAnaHists->GetAnaHist("CMPEtaPhi")->SetAxisLables("#eta","#phi (Rads)");
  _subAnaHists->GetAnaHist("CMUEtaPhi")->SetAxisLables("#eta","#phi (Rads)");
  _subAnaHists->GetAnaHist("CMUPEtaPhi")->SetAxisLables("#eta","#phi (Rads)");
  _subAnaHists->GetAnaHist("dilDphi")->SetAxisLables("Lepton Angular Separation (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiee")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiemu")->SetAxisLables("e #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphimumu")->SetAxisLables("#mu #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("JetNumBeforeCut")->SetAxisLables("Number of Jets","Events");
  _subAnaHists->GetAnaHist("JetNumBeforeCut_ee")->SetAxisLables("Number of Jets","Events");
  _subAnaHists->GetAnaHist("JetNumBeforeCut_emu")->SetAxisLables("Number of Jets","Events");
  _subAnaHists->GetAnaHist("JetNumBeforeCut_mumu")->SetAxisLables("Number of Jets","Events");
  _subAnaHists->GetAnaHist("MissingETMag")->SetAxisLables("Missing Et (GeV)","Events");
  _subAnaHists->GetAnaHist("MissingETMag_ee")->SetAxisLables("Missing Et (GeV)","Events");
  _subAnaHists->GetAnaHist("MissingETMag_emu")->SetAxisLables("Missing Et (GeV)","Events");
  _subAnaHists->GetAnaHist("MissingETMag_mumu")->SetAxisLables("Missing Et (GeV)","Events");
  _subAnaHists->GetAnaHist("dilmassmax")->SetAxisLables("Dilepton Mass (GeV)","Events");
  _subAnaHists->GetAnaHist("dilmassmax_ee")->SetAxisLables("Dilepton Mass (GeV)","Events");
  _subAnaHists->GetAnaHist("dilmassmax_emu")->SetAxisLables("Dilepton Mass (GeV)","Events");
  _subAnaHists->GetAnaHist("dilmassmax_mumu")->SetAxisLables("Dilepton Mass (GeV)","Events");
  _subAnaHists->GetAnaHist("TotalEt")->SetAxisLables("Total Et (GeV)","Events");
  _subAnaHists->GetAnaHist("TotalEt_ee")->SetAxisLables("Total Et (GeV)","Events");
  _subAnaHists->GetAnaHist("TotalEt_emu")->SetAxisLables("Total Et (GeV)","Events");
  _subAnaHists->GetAnaHist("TotalEt_mumu")->SetAxisLables("Total Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_anyN")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_anyN_ee")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_anyN_emu")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_anyN_mumu")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_1Jet")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_1Jet_ee")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_1Jet_emu")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_1Jet_mumu")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_2Jet")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_2Jet_ee")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_2Jet_emu")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("LeadJetEt_2Jet_mumu")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("SecondJetEt_2Jet")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("SecondJetEt_2Jet_ee")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("SecondJetEt_2Jet_emu")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("SecondJetEt_2Jet_mumu")->SetAxisLables("Jet Et (GeV)","Events");
  _subAnaHists->GetAnaHist("TomdilDphi")->SetAxisLables("Lepton #Delta#phi (Rads)","Events");

  _subAnaHists->GetAnaHist("dilDphi_JESup")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiee_JESup")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiemu_JESup")->SetAxisLables("e #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphimumu_JESup")->SetAxisLables("#mu #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphi_JESdn")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiee_JESdn")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiemu_JESdn")->SetAxisLables("e #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphimumu_JESdn")->SetAxisLables("#mu #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphi_LESup")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiee_LESup")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiemu_LESup")->SetAxisLables("e #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphimumu_LESup")->SetAxisLables("#mu #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphi_LESdn")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiee_LESdn")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiemu_LESdn")->SetAxisLables("e #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphimumu_LESdn")->SetAxisLables("#mu #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphi_Fakeup")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiee_Fakeup")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiemu_Fakeup")->SetAxisLables("e #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphimumu_Fakeup")->SetAxisLables("#mu #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphi_Fakedn")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiee_Fakedn")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiemu_Fakedn")->SetAxisLables("e #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphimumu_Fakedn")->SetAxisLables("#mu #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphi_PDFup")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiee_PDFup")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiemu_PDFup")->SetAxisLables("e #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphimumu_PDFup")->SetAxisLables("#mu #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphi_PDFdn")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiee_PDFdn")->SetAxisLables("e e #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphiemu_PDFdn")->SetAxisLables("e #mu #Delta#phi (Rads)","Events");
  _subAnaHists->GetAnaHist("dilDphimumu_PDFdn")->SetAxisLables("#mu #mu #Delta#phi (Rads)","Events");

}

void TDSubAna::GetPDFError(bool high, char * histtofillname, char * weighthist2D, char * midvaluehist)//TD high=true for upper limit
{
  if(_datasource=="data" || _datasource=="fake")
  {
    *((TH1F*)(_subAnaHists->GetAnaHist(histtofillname)->GetHist()))=(*((TH1F*)(_subAnaHists->GetAnaHist(midvaluehist)->GetHist())));
    return;
  }
  TH1 * histtofill=_subAnaHists->GetAnaHist(histtofillname)->GetHist();
  for(int bin=0; bin!=histtofill->GetNbinsX(); ++bin)
  {
    double sumsq=0;
    double max1=0;
    double midvalue=((TH1F*)(_subAnaHists->GetAnaHist(midvaluehist)->GetHist()))->GetBinContent(bin+1);
    if(midvalue==0)
    {
      ((TH1F*)(_subAnaHists->GetAnaHist(histtofillname)->GetHist()))->SetBinContent(bin+1, midvalue);
      continue;
    }
    std::cout<<"midval "<<midvalue<<" bin "<<bin<<std::endl;
    for(int set=0; set<40; set+=2)
    {
      double errorval1=((TH2F*)(_subAnaHists->GetAnaHist(weighthist2D)->GetHist()))->GetBinContent(bin+1, set+1);
      //std::cout<<"errorval1 "<<errorval1<<std::endl;
      double errorval2=((TH2F*)(_subAnaHists->GetAnaHist(weighthist2D)->GetHist()))->GetBinContent(bin+1, set+2);
      //std::cout<<"errorval2 "<<errorval2<<std::endl;
//      if(high) max1=max( (_pdfset[set]-_pdfset[5]), (_pdfset[set+1]-_pdfset[5]) );
//      if(!high)max1=max( (_pdfset[5]-_pdfset[set]), (_pdfset[5]-_pdfset[set+1]) );

      if(high)  max1=max( (errorval1-midvalue)/midvalue, (errorval2-midvalue)/midvalue );//sumsq+=max1*max1;}
      if(!high) max1=min( (errorval1-midvalue)/midvalue, (errorval2-midvalue)/midvalue );//sumsq-=max1*max1;}
      //if(max1>0)
      sumsq+=max1*max1;
    }
    std::cout<<"error: "<<sqrt(sumsq)*midvalue<<std::endl;
    if(high)((TH1F*)(_subAnaHists->GetAnaHist(histtofillname)->GetHist()))->SetBinContent(bin+1, (1+sqrt(sumsq))*midvalue);
    if(!high)((TH1F*)(_subAnaHists->GetAnaHist(histtofillname)->GetHist()))->SetBinContent(bin+1, (1-sqrt(sumsq))*midvalue);
  }
//  return sqrt(sumsq); 
}

double TDSubAna::GetNNCutVal(int cutstrength, int intdiltype)
{
  if(cutstrength==0)
  {
    if(intdiltype==0)return 0.975;
    if(intdiltype==1)return 0.6;
    if(intdiltype==2)return 0.975;
  }else if(cutstrength==1)
  {
    if(intdiltype==0)return 0.985;
    if(intdiltype==1)return 0.8;
    if(intdiltype==2)return 0.985;
  }else if(cutstrength==2)
  {
    if(intdiltype==0)return 0.995;
    if(intdiltype==1)return 0.95;
    if(intdiltype==2)return 0.995;
  }else if(cutstrength==3)
  {
    if(intdiltype==0)return 0.997;
    if(intdiltype==1)return 0.97;
    if(intdiltype==2)return 0.997;
  }else if(cutstrength==4)
  {
    if(intdiltype==0)return 0.9985;
    if(intdiltype==1)return 0.985;
    if(intdiltype==2)return 0.9985;
  }else if(cutstrength==5)
  {
    if(intdiltype==0)return 1.0;
    if(intdiltype==1)return 0.99;
    if(intdiltype==2)return 1.0;
  } 
}

double TDSubAna::GetNNsdyrbCutVal(int intdiltype)
{
  if(_higgsMass==110)
  {
    if(intdiltype==0)return 0.997;
    if(intdiltype==1)return 0.977;
    if(intdiltype==2)return 0.997;
  }else if(_higgsMass==120)
  {
    if(intdiltype==0)return 0.999;
    if(intdiltype==1)return 0.965;
    if(intdiltype==2)return 0.991;
  }else if(_higgsMass==130)
  {
    if(intdiltype==0)return 0.993;
    if(intdiltype==1)return 0.981;
    if(intdiltype==2)return 0.999;
  }else if(_higgsMass==140)
  {
    if(intdiltype==0)return 0.997;
    if(intdiltype==1)return 0.997;
    if(intdiltype==2)return 0.997;
  }else if(_higgsMass==150)
  {
    if(intdiltype==0)return 0.995;
    if(intdiltype==1)return 0.999;
    if(intdiltype==2)return 0.997;
  }else if(_higgsMass==160)
  {
    if(intdiltype==0)return 0.995;
    if(intdiltype==1)return 0.997;
    if(intdiltype==2)return 0.991;
  }else if(_higgsMass==170)
  {
    if(intdiltype==0)return 0.995;
    if(intdiltype==1)return 0.995;
    if(intdiltype==2)return 0.995;
  }else if(_higgsMass==180)
  {
    if(intdiltype==0)return 0.991;
    if(intdiltype==1)return 0.973;
    if(intdiltype==2)return 0.999;
  }else if(_higgsMass==190)
  {
    if(intdiltype==0)return 0.987;
    if(intdiltype==1)return 0.983;
    if(intdiltype==2)return 0.991;
  }else if(_higgsMass==200)
  {
    if(intdiltype==0)return 0.991;
    if(intdiltype==1)return 0.995;
    if(intdiltype==2)return 0.999;
  } 
}

void TDSubAna::SetEpsVtx()
{
 _runint[0]=141543;
 _runint[1]=152945;
 _runint[2]=159602;
 _runint[3]=162100;
 _runint[4]=165063;
 _runint[5]=167783;
 _runint[6]=176568;
 _runint[7]=178759;
 _runint[8]=183697;
 _runint[9]=185010;
 _runint[10]=190697;
 _runint[11]=192401;
 _runint[12]=194101;
 _runint[13]=195801;
 _runint[14]=197501;
 _runint[15]=199201;
 _runint[16]=202301;
 _runint[17]=204549;
 _runint[18]=206169;
 _runint[19]=208792;
 _runint[20]=217990;
 _runint[21]=222529;
 _runint[22]=228664;
 _runint[23]=233133;
 _runint[24]=237845;
 _runint[25]=240788;
 _runint[26]=242055;
 
 _epsvtxvals[0]=0.954;
 _epsvtxvals[1]=0.947;
 _epsvtxvals[2]=0.950;
 _epsvtxvals[3]=0.950;
 _epsvtxvals[4]=0.942;
 _epsvtxvals[5]=0.945;
 _epsvtxvals[6]=0.947;
 _epsvtxvals[7]=0.949;
 _epsvtxvals[8]=0.962;
 _epsvtxvals[9]=0.960;
 _epsvtxvals[10]=0.956;
 _epsvtxvals[11]=0.961;
 _epsvtxvals[12]=0.961;
 _epsvtxvals[13]=0.959;
 _epsvtxvals[14]=0.961;
 _epsvtxvals[15]=0.959;
 _epsvtxvals[16]=0.959;
 _epsvtxvals[17]=0.968;
 _epsvtxvals[18]=0.962;
 _epsvtxvals[19]=0.967;
 _epsvtxvals[20]=0.966;
 _epsvtxvals[21]=0.969;
 _epsvtxvals[22]=0.969;
 _epsvtxvals[23]=0.971;
 _epsvtxvals[24]=0.971;
 _epsvtxvals[25]=0.975;
}

double  TDSubAna::GetEpsVtx(int run)
{
  for(int i=0; i<26; i++)
  {
    if(run>=_runint[i+1]) continue;
    return _epsvtxvals[i];
  }

  std::cout<<" Please update the GetEpsVtx data "<<std::endl;
  return _epsvtxvals[25] ;
}

double TDSubAna::GetNewDYWeight(int leptype, int run)
{
  if(run>0 && run<186599)
  switch(leptype)
  {
    case(0):return 1.001; break;
    case(1):return 1.027;break;
    case(2):return 0;break;
    case(3):return 1.004;break;
    case(4):return 1.012;break;
    case(5):return 0.936;break;
    case(6):return 0.974;break;
    case(7):return 0.706;break;
    case(8):return 1;break;
    case(9):return 1;break;
    case(10):return 1;break;
    case(11):return 1;break;
    case(12):
	    default:
	     return 1;break;
  }
  else if(run<203799)
    switch(leptype) 
    { 
      case(0):return 1.004; break;
      case(1):return 0.990;break;
      case(2):return 0;break;
      case(3):return 0.997;break;
      case(4):return 1.011;break;
      case(5):return 0.991;break;
      case(6):return 0.990;break;
      case(7):return 0.988;break;
      case(8):return 1;break;
      case(9):return 1;break;
      case(10):return 1;break;
      case(11):return 1;break;
      case(12):
      default:
	       return 1;break;
    }
  else
    switch(leptype) 
    { 
      case(0):return 0.997; break;
      case(1):return 0.986;break;
      case(2):return 0;break;
      case(3):return 0.999;break;
      case(4):return 1.013;break;
      case(5):return 0.982;break;
      case(6):return 0.973;break;
      case(7):return 0.978;break;
      case(8):return 1;break;
      case(9):return 1;break;
      case(10):return 1;break;
      case(11):return 1;break;
      case(12):
      default:
	       return 1;break;
    }
}