#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 = new TMVA::Reader( *inputVars );
       reader2 = new TMVA::Reader( *inputVars );

       string dir    = "./weights/";
       char prefix[200];
       sprintf(prefix,"TMV1_%d", (int)_higgsMass);
       //string prefix = "TMVAnalysis";
       //string prefix =buffer;

   
   if (Use_Cuts)
    reader->BookMVA( "TMVA::Types::Cuts",         dir + prefix + "_Cuts.weights"     );
  if (Use_Likelihood)
    reader->BookMVA( "TMVA::Types::Likelihood",   dir + prefix + "_Likelihood.weights.root"  );
  if (Use_LikelihoodD)
    reader->BookMVA( "TMVA::Types::LikelihoodD",   dir + prefix + "_LikelihoodD.weights.root" );
  if (Use_PDERS)
    reader->BookMVA( "TMVA::Types::PDERS",        dir + prefix + "_PDERS.weights.root"  );
  if (Use_HMatrix)
    reader->BookMVA( "TMVA::Types::HMatrix",      dir + prefix + "_HMatrix.weights"  );
  if (Use_Fisher)
    reader->BookMVA( "TMVA::Types::Fisher",       dir + prefix + "_Fisher.weights"   );
  if (Use_CFMlpANN)
    reader->BookMVA( "TMVA::Types::CFMlpANN",     dir + prefix + "_CFMlpANN.weights" );
  if (Use_TMlpANN)
    reader->BookMVA( "TMVA::Types::TMlpANN",      dir + prefix + "_TMlpANNBFGS.weights.txt");//"_TMlpANN.weights"  );//quick fix for now
  if (Use_BDT)
    reader->BookMVA( "TMVA::Types::BDT",          dir + prefix + "_BDTGini.weights"  );

  sprintf(prefix,"TMV2_%d", (int)_higgsMass);
  //string prefix =buffer;

  if (Use_Cuts)
    reader2->BookMVA( "TMVA::Types::Cuts",         dir + prefix + "_Cuts.weights"     );
  if (Use_Likelihood)
    reader2->BookMVA( "TMVA::Types::Likelihood",   dir + prefix + "_Likelihood.weights.root"  );
  if (Use_LikelihoodD)
    reader2->BookMVA( "TMVA::Types::LikelihoodD",   dir + prefix + "_LikelihoodD.weights.root" );
  if (Use_PDERS)
    reader2->BookMVA( "TMVA::Types::PDERS",        dir + prefix + "_PDERS.weights.root"  );
  if (Use_HMatrix)
    reader2->BookMVA( "TMVA::Types::HMatrix",      dir + prefix + "_HMatrix.weights"  );
  if (Use_Fisher)
    reader2->BookMVA( "TMVA::Types::Fisher",       dir + prefix + "_Fisher.weights"   );
  if (Use_CFMlpANN)
    reader2->BookMVA( "TMVA::Types::CFMlpANN",     dir + prefix + "_CFMlpANN.weights" );
  if (Use_TMlpANN)
    reader2->BookMVA( "TMVA::Types::TMlpANN",      dir + prefix + "_TMlpANNBFGS.weights.txt");//"_TMlpANN.weights"  );
  if (Use_BDT)
    reader2->BookMVA( "TMVA::Types::BDT",          dir + prefix + "_BDTGini.weights"  );
  initTMVA=true;
}

bool TDSubAna::evaluateTMVA(double weight, TDSubAnaCutInterface * cutInterface)
{
  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);
  //varValues.push_back(_NNdiltype);

  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 (Use_Likelihood )  MVA_Likelihood  = reader->EvaluateMVA( varValues, "TMVA::Types::Likelihood" );
 if (Use_LikelihoodD)  MVA_LikelihoodD = reader->EvaluateMVA( varValues, "TMVA::Types::LikelihoodD");//"1                       );
 if (Use_PDERS      )  MVA_PDERS       = reader->EvaluateMVA( varValues, "TMVA::Types::PDERS"      );
 if (Use_HMatrix    )  MVA_HMatrix     = reader->EvaluateMVA( varValues, "TMVA::Types::HMatrix"    );
 if (Use_Fisher     )  MVA_Fisher      = reader->EvaluateMVA( varValues, "TMVA::Types::Fisher"     );
 if (Use_CFMlpANN   )  MVA_CFMlpANN    = reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   );
 if (Use_TMlpANN    )  MVA_TMlpANN     = reader->EvaluateMVA( varValues, "TMVA::Types::TMlpANN"    );
 if (Use_BDT        )  MVA_BDT         = reader->EvaluateMVA( varValues, "TMVA::Types::BDT"        );

 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 (Use_Likelihood )  MVA_Likelihood2  = reader2->EvaluateMVA( varValues, "TMVA::Types::Likelihood" );//TMVA::Types::Likelihood );
 if (Use_LikelihoodD)  MVA_LikelihoodD2 = reader2->EvaluateMVA( varValues, "TMVA::Types::LikelihoodD");//1                       );
 if (Use_PDERS      )  MVA_PDERS2       = reader2->EvaluateMVA( varValues, "TMVA::Types::PDERS"      );//TMVA::Types::PDERS      );
 if (Use_HMatrix    )  MVA_HMatrix2     = reader2->EvaluateMVA( varValues, "TMVA::Types::HMatrix"    );//TMVA::Types::HMatrix    );
 if (Use_Fisher     )  MVA_Fisher2      = reader2->EvaluateMVA( varValues, "TMVA::Types::Fisher"     );//TMVA::Types::Fisher     );
 if (Use_CFMlpANN   )  MVA_CFMlpANN2    = reader2->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   );//TMVA::Types::CFMlpANN   );
 if (Use_TMlpANN    )  MVA_TMlpANN2     = reader2->EvaluateMVA( varValues, "TMVA::Types::TMlpANN"    );//TMVA::Types::TMlpANN    );
 if (Use_BDT        )  MVA_BDT2         = reader2->EvaluateMVA( varValues, "TMVA::Types::BDT"        );//TMVA::Types::BDT        );

  if (Use_Likelihood ) HistFill("MVA_Likelihood" , MVA_Likelihood , weight);
  if (Use_LikelihoodD) HistFill("MVA_LikelihoodD", MVA_LikelihoodD, weight);
  if (Use_PDERS      ) HistFill("MVA_PDERS"      , MVA_PDERS      , weight);
  if (Use_HMatrix    ) HistFill("MVA_HMatrix"    , MVA_HMatrix    , weight);
  if (Use_Fisher     ) HistFill("MVA_Fisher"     , MVA_Fisher     , weight);
  if (Use_CFMlpANN   ) HistFill("MVA_CFMlpANN"   , MVA_CFMlpANN   , weight);
  if (Use_TMlpANN    ) HistFill("MVA_TMlpANN"    , MVA_TMlpANN    , weight);
  if (Use_BDT        ) HistFill("MVA_BDT"        , MVA_BDT        , weight);

  if (Use_Likelihood ) HistFill("MVA2_Likelihood" , MVA_Likelihood2 , weight);
  if (Use_LikelihoodD) HistFill("MVA2_LikelihoodD", MVA_LikelihoodD2, weight);
  if (Use_PDERS      ) HistFill("MVA2_PDERS"      , MVA_PDERS2      , weight);
  if (Use_HMatrix    ) HistFill("MVA2_HMatrix"    , MVA_HMatrix2    , weight);
  if (Use_Fisher     ) HistFill("MVA2_Fisher"     , MVA_Fisher2     , weight);
  if (Use_CFMlpANN   ) HistFill("MVA2_CFMlpANN"   , MVA_CFMlpANN2   , weight);
  if (Use_TMlpANN    ) HistFill("MVA2_TMlpANN"    , MVA_TMlpANN2    , weight);
  if (Use_BDT        ) HistFill("MVA2_BDT"        , MVA_BDT2        , weight);

  if(Use_CFMlpANN    && reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   )<0.002)HistFill("MVA2_CFMlpANN_cut2_0"   , MVA_CFMlpANN2   , weight);
  if(Use_CFMlpANN    && reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   )<0.004)HistFill("MVA2_CFMlpANN_cut2_1"   , MVA_CFMlpANN2   , weight);
  if(Use_CFMlpANN    && reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   )<0.006)HistFill("MVA2_CFMlpANN_cut2_2"   , MVA_CFMlpANN2   , weight);
  if(Use_CFMlpANN    && reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   )<0.008)HistFill("MVA2_CFMlpANN_cut2_3"   , MVA_CFMlpANN2   , weight);
  if(Use_CFMlpANN    && reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   )<0.01)HistFill("MVA2_CFMlpANN_cut2_4"   , MVA_CFMlpANN2   , weight);
  if(Use_CFMlpANN    && reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   )<0.012)HistFill("MVA2_CFMlpANN_cut2_5"   , MVA_CFMlpANN2   , weight);
  if(Use_CFMlpANN    && reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   )<0.014)HistFill("MVA2_CFMlpANN_cut2_6"   , MVA_CFMlpANN2   , weight);
  if(Use_CFMlpANN    && reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   )<0.016)HistFill("MVA2_CFMlpANN_cut2_7"   , MVA_CFMlpANN2   , weight);
  if(Use_CFMlpANN    && reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   )<0.018)HistFill("MVA2_CFMlpANN_cut2_8"   , MVA_CFMlpANN2   , weight);
  if(Use_CFMlpANN    && reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   )<0.02)HistFill("MVA2_CFMlpANN_cut2_9"   , MVA_CFMlpANN2   , weight);

  /*if(Use_HMatrix     && reader->EvaluateMVA( varValues, "TMVA::Types::HMatrix"    )<0.5) 
  if(Use_BDT         && reader->EvaluateMVA( varValues, "TMVA::Types::BDT"        )<0.5)
  if(Use_PDERS       && reader->EvaluateMVA( varValues, "TMVA::Types::PDERS"      )<0.1)*/
  
  if(//(Use_Likelihood  && reader->EvaluateMVA( varValues, "TMVA::Types::Likelihood" )<0.2) &&
     //(Use_LikelihoodD && reader->EvaluateMVA( varValues, "TMVA::Types::LikelihoodD" )<0.4) &&
     (Use_HMatrix     && reader->EvaluateMVA( varValues, "TMVA::Types::HMatrix"    )<0.5) &&
     //(Use_Fisher      && reader->EvaluateMVA( varValues, "TMVA::Types::Fisher"     )<1)&&
     (Use_BDT         && reader->EvaluateMVA( varValues, "TMVA::Types::BDT"        )<0.5)&&
     (Use_PDERS       && reader->EvaluateMVA( varValues, "TMVA::Types::PDERS"      )<0.1)&&
     (Use_CFMlpANN    && reader->EvaluateMVA( varValues, "TMVA::Types::CFMlpANN"   )<0.1))
  {
    if (Use_Likelihood ) HistFill("MVA_Likelihood_cut1" , MVA_Likelihood , weight);
    if (Use_LikelihoodD) HistFill("MVA_LikelihoodD_cut1", MVA_LikelihoodD, weight);
    if (Use_PDERS      ) HistFill("MVA_PDERS_cut1"      , MVA_PDERS      , weight);
    if (Use_HMatrix    ) HistFill("MVA_HMatrix_cut1"    , MVA_HMatrix    , weight);
    if (Use_Fisher     ) HistFill("MVA_Fisher_cut1"     , MVA_Fisher     , weight);
    if (Use_CFMlpANN   ) HistFill("MVA_CFMlpANN_cut1"   , MVA_CFMlpANN   , weight);
    if (Use_TMlpANN    ) HistFill("MVA_TMlpANN_cut1"    , MVA_TMlpANN    , weight);
    if (Use_BDT        ) HistFill("MVA_BDT_cut1"        , MVA_BDT        , weight);

    if (Use_Likelihood ) HistFill("MVA2_Likelihood_cut1" , MVA_Likelihood2 , weight);
    if (Use_LikelihoodD) HistFill("MVA2_LikelihoodD_cut1", MVA_LikelihoodD2, weight);
    if (Use_PDERS      ) HistFill("MVA2_PDERS_cut1"      , MVA_PDERS2      , weight);
    if (Use_HMatrix    ) HistFill("MVA2_HMatrix_cut1"    , MVA_HMatrix2    , weight);
    if (Use_Fisher     ) HistFill("MVA2_Fisher_cut1"     , MVA_Fisher2     , weight);
    if (Use_CFMlpANN   ) HistFill("MVA2_CFMlpANN_cut1"   , MVA_CFMlpANN2   , weight);
    if (Use_TMlpANN    ) HistFill("MVA2_TMlpANN_cut1"    , MVA_TMlpANN2    , weight);
    if (Use_BDT        ) HistFill("MVA2_BDT_cut1"        , MVA_BDT2        , weight);

    if((Use_Likelihood  && reader2->EvaluateMVA( varValues, TMVA::Types::Likelihood )<0.4) &&
     (Use_LikelihoodD && reader2->EvaluateMVA( varValues, 1                         )<0.6) &&
//     (Use_HMatrix     && reader2->EvaluateMVA( varValues, TMVA::Types::HMatrix    )<0.5) &&
     (Use_Fisher      && reader2->EvaluateMVA( varValues, TMVA::Types::Fisher       )<1  ) &&
     (Use_PDERS       &&  reader2->EvaluateMVA( varValues, TMVA::Types::PDERS       )<0.5  ) &&
//&&     (Use_CFMlpANN    && reader2->EvaluateMVA( varValues, TMVA::Types::CFMlpANN   )<0.7))
     (Use_Fisher      && reader2->EvaluateMVA( varValues, TMVA::Types::BDT          )<0.6  ))
    {
      if (Use_Likelihood ) HistFill("MVA2_Likelihood_cut2" , MVA_Likelihood2 , weight);
      if (Use_LikelihoodD) HistFill("MVA2_LikelihoodD_cut2", MVA_LikelihoodD2, weight);
      if (Use_PDERS      ) HistFill("MVA2_PDERS_cut2"      , MVA_PDERS2      , weight);
      if (Use_HMatrix    ) HistFill("MVA2_HMatrix_cut2"    , MVA_HMatrix2    , weight);
      if (Use_Fisher     ) HistFill("MVA2_Fisher_cut2"     , MVA_Fisher2     , weight);
      if (Use_CFMlpANN   ) HistFill("MVA2_CFMlpANN_cut2"   , MVA_CFMlpANN2   , weight);
      if (Use_TMlpANN    ) HistFill("MVA2_TMlpANN_cut2"    , MVA_TMlpANN2    , weight);
      if (Use_BDT        ) HistFill("MVA2_BDT_cut2"        , MVA_BDT2        , weight);
      if (MVA_CFMlpANN2<0.5)HistFill("dilDphi_TMVA", dildphi, weight);
    }
    
    //_descVar2->Fill();
//    evaluateTMVA2(weight, cutInterface);
    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);
}

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;
  }
//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
  _fakevectord->push_back((TH1F*)fRootFile.Get("CEM_d")->Clone());
  _fakevectord->push_back((TH1F*)fRootFile.Get("PHX_d")->Clone());
  _fakevectord->push_back(NULL);//TD PEM - should never be called!!!
  _fakevectord->push_back((TH1F*)fRootFile.Get("CMUP_d")->Clone());
  _fakevectord->push_back((TH1F*)fRootFile.Get("CMX_d")->Clone());
  _fakevectord->push_back((TH1F*)fRootFile.Get("CMU_d")->Clone());
  _fakevectord->push_back((TH1F*)fRootFile.Get("CMP_d")->Clone());
  _fakevectord->push_back((TH1F*)fRootFile.Get("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("CEM_h")->Clone());
  _fakevectorh->push_back((TH1F*)fRootFile.Get("PHX_h")->Clone());
  _fakevectorh->push_back(NULL);//TD PEM - shoulh never be calleh!!!
  _fakevectorh->push_back((TH1F*)fRootFile.Get("CMUP_h")->Clone());
  _fakevectorh->push_back((TH1F*)fRootFile.Get("CMX_h")->Clone());
  _fakevectorh->push_back((TH1F*)fRootFile.Get("CMU_h")->Clone());
  _fakevectorh->push_back((TH1F*)fRootFile.Get("CMP_h")->Clone());
  _fakevectorh->push_back((TH1F*)fRootFile.Get("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("CEM_i")->Clone());
  _fakevectori->push_back((TH1F*)fRootFile.Get("PHX_i")->Clone());
  _fakevectori->push_back(NULL);//TD PEM - shouli never be callei!!!
  _fakevectori->push_back((TH1F*)fRootFile.Get("CMUP_i")->Clone());
  _fakevectori->push_back((TH1F*)fRootFile.Get("CMX_i")->Clone());
  _fakevectori->push_back((TH1F*)fRootFile.Get("CMU_i")->Clone());
  _fakevectori->push_back((TH1F*)fRootFile.Get("CMP_i")->Clone());
  _fakevectori->push_back((TH1F*)fRootFile.Get("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
    {
  //    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();
}