#define opt_disc_cxx
#include "opt_disc.h"
#include <TROOT.h>
#include <TChain.h>
#include <TFile.h>
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TEventList.h"
#include <fstream>
#include <vector>
#include "TMath.h"
#include "TPave.h"
#include <map>
#include <vector>
#include <iostream>
#include "TString.h"

using namespace std;

Double_t        P_tb;
Double_t        P_tq;
Double_t        P_tbg;
Double_t        P_tqb;
Double_t        P_tqg;
Double_t        P_wbb;
Double_t        P_wcg;
Double_t        P_wgg;
Double_t        P_wbbg;
Double_t        P_wugg;
Double_t        P_WW;
Double_t        P_WZ;
Double_t        P_lepjets_2jet;
Double_t        P_ggg;
Double_t        P_lepjets_3jet_1;
Double_t        P_lepjets_3jet_2;
Double_t        c1;
Double_t        c2;
Double_t        c3;
Double_t        c4;
Double_t        c5;
Double_t        c6;
Double_t        c7;
Double_t        acc_s;
Double_t        acc_t;
TFile*          outfile;
TH1F *hist_st;
TH1F *hist_btag_st;
TH1F *hist_st1; 
TH1F *hist_st2;
TH1F *hist_st3;
TH1F *hist_st4;
TH1F *hist_st5;
TH1F *hist_st6;
TH1F *hist_s;
TH1F *hist_s1;
TH1F *hist_s2;
TH1F *hist_s3;
TH1F *hist_s4;
TH1F *hist_s5;
TH1F *hist_s6;
TH1F *hist_t;
TH1F *hist_t1;
TH1F *hist_t2;
TH1F *hist_t3;
TH1F *hist_t4;
TH1F *hist_t5;
TH1F *hist_t6;
TH1F *hist_r;
TH1F *hist_r_2d;
TH2F *hist_st_2d;
TH2F *hist_bst_2d;
TH1F *hist_bst2;
TH1F *hist_bs2;
TH1F *hist_bt2;
TTree *TopologicalVariables_new;

Int_t njet;

double _evtwgt;
double _Prob_tb;
double _Prob_tbg;
double _Prob_tq;
double _Prob_tqb;
double _Prob_tqg;
double _Prob_wugg;
double _Prob_wbb;
double _Prob_wbbg;
double _Prob_wcg;
double _Prob_wgg;
double _Prob_WW;
double _Prob_WZ;
double _Prob_lepjets_2jet;
double _Prob_ggg;
double _Prob_lepjets_3jet_1;
double _Prob_lepjets_3jet_2;
double _LeadingJetBTagNN;
double _SecondJetBTagNN;
double _Jet1BTagNN;
double _Jet2BTagNN;
double _disc_tb_btag;
double _disc_tq_btag;
double _disc_tb_tq_btag;

double bwgt;
double not_bwgt;
double evtwgt_scale;

void opt_disc::Init(string inputrootfile, string outputrootfilename, int numjets, double y_tb, double y_tqb, double y_wbb, double y_wcc, double y_wlp, double y_qcd, double y_diboson, double y_ttbar, double scale) {

  cout << " file " << inputrootfile.c_str() << endl;

  TFile *MEFile = new TFile(inputrootfile.c_str());

  TTree *METree = (TTree*)MEFile->Get("TopologicalVariables");
  Init(METree);

  double tot=y_wbb+y_wcc+y_wlp+y_qcd+y_diboson+y_ttbar;
  double tot_sig=y_tb+y_tqb;
  acc_s=y_tb/tot_sig;
  acc_t=y_tqb/tot_sig;

 
  c1=y_wbb/tot+.5*y_wcc/tot;
  c2=.5*y_wcc/tot;
  c3=y_wlp/tot+.25*y_qcd/tot;
  c4=.5*y_diboson/tot;
  c5=.5*y_diboson/tot;
  c6=.75*y_qcd/tot;
  c7=y_ttbar/tot;

  evtwgt_scale=scale;

  cout << acc_s << " " << acc_t << " " << c1 << " " << c2 << " " << c3 << " " << c4 << " " << c5 << " " << c6 << " " << c7 << endl; 

  njet=numjets;

  outfile = new TFile(outputrootfilename.c_str(), "RECREATE");
  hist_st = new TH1F( "DiscriminantST_all", "DiscriminantS_all", 500, 0,1);
  hist_st1 = new TH1F( "DiscriminantST_all_1000", "DiscriminantST_all_1000", 1000, 0,1);
  hist_st2 = new TH1F( "DiscriminantST_all_100", "DiscriminantST_all_100", 100, 0,1);
  hist_st3 = new TH1F( "DiscriminantST_all_10", "DiscriminantST_all_10", 10, 0,1);
  hist_st4 = new TH1F( "DiscriminantST_all_20", "DiscriminantST_all_20", 20, 0,1);
  hist_st5 = new TH1F( "DiscriminantST_all_50", "DiscriminantST_all_50", 50, 0,1);
  hist_st6 = new TH1F( "DiscriminantST_all_25", "DiscriminantST_all_25", 25, 0,1);
  hist_s = new TH1F( "DiscriminantS_all", "DiscriminantS_all", 500, 0,1);
  hist_s1 = new TH1F( "DiscriminantS_all_1000", "DiscriminantS_all_1000", 1000, 0,1);
  hist_s2 = new TH1F( "DiscriminantS_all_100", "DiscriminantS_all_100", 100, 0,1);
  hist_s3 = new TH1F( "DiscriminantS_all_10", "DiscriminantS_all_10", 10, 0,1);
  hist_s4 = new TH1F( "DiscriminantS_all_20", "DiscriminantS_all_20", 20, 0,1);
  hist_s5 = new TH1F( "DiscriminantS_all_50", "DiscriminantS_all_50", 50, 0,1);
  hist_s6 = new TH1F( "DiscriminantS_all_25", "DiscriminantS_all_25", 25, 0,1);
  hist_t = new TH1F( "DiscriminantT_all", "DiscriminantT_all", 500, 0,1);
  hist_t1 = new TH1F( "DiscriminantT_all_1000", "DiscriminantT_all_1000", 1000, 0,1);
  hist_t2 = new TH1F( "DiscriminantT_all_100", "DiscriminantT_all_100", 100, 0,1);
  hist_t3 = new TH1F( "DiscriminantT_all_10", "DiscriminantT_all_10", 10, 0,1);
  hist_t4 = new TH1F( "DiscriminantT_all_20", "DiscriminantT_all_20", 20, 0,1);
  hist_t5 = new TH1F( "DiscriminantT_all_50", "DiscriminantT_all_50", 50, 0,1);
  hist_t6 = new TH1F( "DiscriminantT_all_25", "DiscriminantT_all_25", 25, 0,1);
  hist_bst2 = new TH1F( "DiscriminantST_btag_100", "DiscriminantST_btag_100", 100, 0,1);
  hist_bs2 = new TH1F( "DiscriminantS_btag_100", "DiscriminantS_btag_100", 100, 0,1);
  hist_bt2 = new TH1F( "DiscriminantT_btag_100", "DiscriminantT_btag_100", 100, 0,1);


  hist_r = new TH1F( "Discriminant_r_100", "Discriminant_r_100", 100, 0,1);
  hist_r_2d = new TH1F( "Discriminant_r_2d_100", "Discriminant_r_2d_100", 100, 0,1);
  hist_st_2d = new TH2F( "Discriminant_st_2d_100", "Discriminant_st_2d_100", 100, 0,1,100,0,1);
  hist_bst_2d = new TH2F( "Discriminant_btag_st_2d_100", "Discriminant_btag_st_2d_100", 100, 0,1,100,0,1);
  
  METree->SetBranchStatus("Disc_tb_btag", 0);
  METree->SetBranchStatus("Disc_tq_btag", 0);
  METree->SetBranchStatus("Disc_tb_tq_btag", 0);

  //TopologicalVariables_new = new TTree("TopologicalVariables","An example of a ROOT tree");
  TopologicalVariables_new = METree->CloneTree(0);
  /*
  METree->SetBranchStatus("Disc_tb_btag_btag", 1);
  METree->SetBranchStatus("Disc_tq_btag_btag", 1);
  METree->SetBranchStatus("Disc_tb_tq_btag_btag", 1);
  */

  TopologicalVariables_new->Branch("Disc_tb_btag", &_disc_tb_btag, "Disc_tb_btag/D");
  TopologicalVariables_new->Branch("Disc_tb_tq_btag", &_disc_tb_tq_btag, "Disc_tb_tq_btag/D");
  TopologicalVariables_new->Branch("Disc_tq_btag", &_disc_tq_btag, "Disc_tq_btag/D");

  /*
 TopologicalVariables_new = new TTree("TopologicalVariables","An example of a ROOT tree");

  TopologicalVariables_new->Branch("EventWeight", &_evtwgt, "EventWeight/D");
  TopologicalVariables_new->Branch("Probability_tb_btag", &_Prob_tb, "Probability_tb_btag/D");
  TopologicalVariables_new->Branch("Probability_tq_btag", &_Prob_tq, "Probability_tq_btag/D");
  TopologicalVariables_new->Branch("Probability_tbg_btag", &_Prob_tbg, "Probability_tbg_btag/D");
  TopologicalVariables_new->Branch("Probability_tqg_btag", &_Prob_tqg, "Probability_tqg_btag/D");
  TopologicalVariables_new->Branch("Probability_tqb_btag", &_Prob_tqb, "Probability_tqb_btag/D");
  TopologicalVariables_new->Branch("Probability_wbb_btag", &_Prob_wbb, "Probability_wbb_btag/D");
  TopologicalVariables_new->Branch("Probability_wcg_btag", &_Prob_wcg, "Probability_wcg_btag/D");
  TopologicalVariables_new->Branch("Probability_wgg_btag", &_Prob_wgg, "Probability_wgg_btag/D");
  TopologicalVariables_new->Branch("Probability_wbbg_btag", &_Prob_wbbg, "Probability_wbbg_btag/D");
  TopologicalVariables_new->Branch("Probability_WW_btag", &_Prob_WW, "Probability_WW_btag/D");
  TopologicalVariables_new->Branch("Probability_WZ_btag", &_Prob_WZ, "Probability_WZ_btag/D");
  TopologicalVariables_new->Branch("Probability_lepjets_2jet_btag", &_Prob_lepjets_2jet, "Probability_lepjets_2jet_btag/D");
  TopologicalVariables_new->Branch("Probability_ggg_btag", &_Prob_ggg, "Probability_ggg_btag/D");
  TopologicalVariables_new->Branch("Probability_lepjets_3jet_1_btag", &_Prob_lepjets_3jet_1, "Probability_lepjets_3jet_1_btag/D");
  TopologicalVariables_new->Branch("Probability_lepjets_3jet_2_btag", &_Prob_lepjets_3jet_2, "Probability_lepjets_3jet_2_btag/D");
  TopologicalVariables_new->Branch("Probability_wugg_btag", &_Prob_wugg, "Probability_wugg_btag/D");
  TopologicalVariables_new->Branch("LeadingBTaggedJetBTagNN", &_LeadingJetBTagNN, "LeadingJetBTagNN/D");
  TopologicalVariables_new->Branch("SecondBTaggedJetBTagNN", &_SecondJetBTagNN, "SecondJetBTagNN/D");
  TopologicalVariables_new->Branch("Jet1BTagNN", &_Jet1BTagNN, "Jet1BTagNN/D");
  TopologicalVariables_new->Branch("Jet2BTagNN", &_Jet2BTagNN, "Jet2BTagNN/D");
 */

  
}




void opt_disc::Loop()
{
//   In a ROOT session, you can do:
//      Root > .L opt_disc.C
//      Root > opt_disc t
//      Root > t.GetEntry(12); // Fill t data members with entry number 12
//      Root > t.Show();       // Show values of entry 12
//      Root > t.Show(16);     // Read and show values of entry 16
//      Root > t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();
   int count=0;
   Int_t nbytes = 0, nb = 0;
   double eventweight;

   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      eventweight = EventWeight*evtwgt_scale;
      P_tb=		      Probability_tb_btag;		 
      P_tq=		      Probability_tq_btag;		 
      P_tbg=		      Probability_tbg_btag;		 
      P_tqb=		      Probability_tqb_btag;		 
      P_tqg=		      Probability_tqg_btag;		 
      P_wbb=		      Probability_wbb_btag;		 
      P_wcg=		      Probability_wcg_btag;		 
      P_wgg=		      Probability_wgg_btag;		 
      P_wbbg=		      Probability_wbbg_btag;		 
      P_wugg=		      Probability_wugg_btag;		 
      P_WW=		      Probability_WW_btag;		 
      P_WZ=		      Probability_WZ_btag;		 
      P_lepjets_2jet=	      Probability_lepjets_2jet_btag;	 
      P_ggg=		      Probability_ggg_btag;		 
      P_lepjets_3jet_1=     Probability_lepjets_3jet_1_btag;
      P_lepjets_3jet_2=     Probability_lepjets_3jet_2_btag;

      /*
      cout << P_tb << " " << P_tq  << " " << P_tbg  << " " << P_tqb  << " ";
      cout <<P_tqg  << " " <<P_wbb  << " " <<P_wcg  << " " <<P_wgg  << " ";
      cout <<P_wbbg  << " "  << " " <<P_wugg  << " " << P_WW  << " ";
      cout <<P_WZ  << " " << P_lepjets_2jet  << " " << P_ggg  << " ";
      cout  << " " << P_lepjets_3jet_1  << " " <<P_lepjets_3jet_2 << endl;
      */

      if(P_tb < 0) P_tb =0.0;
      if(P_tbg < 0) P_tbg =0.0;
      if(P_tq < 0) P_tq =0.0;
      if(P_tqb < 0) P_tqb =0.0;
      if(P_wbb < 0) P_wbb =0.0;
      if(P_wcg < 0) P_wcg =0.0;
      if(P_wgg < 0) P_wgg =0.0;
      if(P_wbbg < 0) P_wbbg =0.0;
      if(P_WW < 0) P_WW =0.0;
      if(P_WZ < 0) P_WZ =0.0;
      if(P_ggg < 0) P_ggg =0.0;
      if(P_lepjets_2jet < 0) P_lepjets_2jet =0.0;
      if(P_lepjets_3jet_1 < 0) P_lepjets_3jet_1 =0.0;
      if(P_lepjets_3jet_2 < 0) P_lepjets_3jet_2 =0.0;


      _evtwgt=EventWeight;
      _Prob_tb =P_tb;
      _Prob_tbg=P_tbg;
      _Prob_tq =P_tq;
      _Prob_tqb=P_tqb;
      _Prob_tqg=P_tqg;
      _Prob_wugg=P_wugg;
      _Prob_wbb=P_wbb;
      _Prob_wbbg=P_wbbg;
      _Prob_wcg=P_wcg;
      _Prob_wgg=P_wgg;
      _Prob_WW=P_WW;
      _Prob_WZ=P_WZ;
      _Prob_lepjets_2jet=P_lepjets_2jet;
      _Prob_ggg=P_ggg;
      _Prob_lepjets_3jet_1=P_lepjets_3jet_1;
      _Prob_lepjets_3jet_2=P_lepjets_3jet_2;
      _LeadingJetBTagNN= LeadingBTaggedJetBTagNN;
      _SecondJetBTagNN=  SecondBTaggedJetBTagNN;
      _Jet1BTagNN=Jet1BTagNN;
      _Jet2BTagNN=Jet2BTagNN;



      double Dst_1 = CalculateD_st(njet);
      double Ds_1 =  CalculateD_s(njet);
      double Dt_1 =  CalculateD_t(njet);
      double r=sqrt(Dst_1*Dst_1+Ds_1*Ds_1+Dt_1*Dt_1);
      double r2d=sqrt(acc_s*acc_s*Ds_1*Ds_1+acc_t*acc_t*Dt_1*Dt_1);
      if (count==0) cout << njet<< " " <<acc_s<< " " <<acc_t<< " " <<c1<< " " <<c2<< " " <<c3<< " " <<c4<< " " <<c5<< " " <<c6<< " " <<c7 << " " << Ds_1 << " " << Dt_1 << " " << Dst_1 << " " << r << endl;
      
      if ((Ds_1 == -1) || (Ds_1 < 0) || (Dt_1 == -1) || (Dt_1 < 0) ) {
      cout << Probability_tb_btag << " " << Probability_tq_btag  << " " << Probability_tbg_btag  << " " << Probability_tqb_btag  << " ";
      cout <<Probability_tqg_btag  << " " <<Probability_wbb_btag  << " " <<Probability_wcg_btag  << " " <<Probability_wgg_btag  << " ";
      cout <<Probability_wbbg_btag  << " "  << " " <<Probability_wugg_btag  << " " << Probability_WW_btag  << " ";
      cout <<Probability_WZ_btag  << " " << Probability_lepjets_2jet_btag  << " " << Probability_ggg_btag  << " ";
      cout  << " " << Probability_lepjets_3jet_1_btag  << " " <<Probability_lepjets_3jet_2_btag << endl;
      string keys = "";
      keys=MCkey;
      cout << " MC Key " << keys << endl;
      }


      int nbtag=0;
      if ((Jet1BTagNN>0.05) && (Jet2BTagNN>0.05)) {
	nbtag=2;
	bwgt= ((Jet1BTagNN-0.8)/0.2)*((Jet2BTagNN-0.8)/0.2);
	not_bwgt=(1-((Jet1BTagNN-0.8)/0.2))*(1-((Jet2BTagNN-0.8)/0.2));
      } else if ((Jet1BTagNN >0.05) && (Jet2BTagNN<0.05)){
	nbtag=1;
	bwgt= ((Jet1BTagNN-0.8)/0.2);
	not_bwgt=(1-bwgt);
      } else if  ((Jet1BTagNN <0.05) & (Jet2BTagNN>0.05)) {
	nbtag=1;
	bwgt= ((Jet2BTagNN-0.8)/0.2);
	not_bwgt=(1-bwgt);
      }

      double Dst_btag_1 = CalculateD_btag_st(njet);
      double Ds_btag_1 =  CalculateD_btag_s(njet);
      double Dt_btag_1 =  CalculateD_btag_t(njet);


      _disc_tb_btag=Ds_1;
      _disc_tq_btag=Dt_1;
      _disc_tb_tq_btag=Dst_1;

      hist_r_2d->Fill(r2d,eventweight);
      hist_r->Fill(r,eventweight);
      hist_s->Fill(Ds_1,eventweight);
      hist_s1->Fill(Ds_1,eventweight);
      hist_s2->Fill(Ds_1,eventweight);
      hist_s3->Fill(Ds_1,eventweight);
      hist_s4->Fill(Ds_1,eventweight);
      hist_s5->Fill(Ds_1,eventweight);
      hist_s6->Fill(Ds_1,eventweight);
      hist_t->Fill(Dt_1,eventweight);
      hist_t1->Fill(Dt_1,eventweight);
      hist_t2->Fill(Dt_1,eventweight);
      hist_t3->Fill(Dt_1,eventweight);
      hist_t4->Fill(Dt_1,eventweight);
      hist_t5->Fill(Dt_1,eventweight);
      hist_t6->Fill(Dt_1,eventweight);
      hist_st->Fill(Dst_1,eventweight);
      hist_st1->Fill(Dst_1,eventweight);
      hist_st2->Fill(Dst_1,eventweight);
      hist_st3->Fill(Dst_1,eventweight);
      hist_st4->Fill(Dst_1,eventweight);
      hist_st5->Fill(Dst_1,eventweight);
      hist_st6->Fill(Dst_1,eventweight);
      hist_st_2d->Fill(Ds_1,Dt_1,eventweight);


      hist_bst2->Fill(Dst_btag_1,eventweight);
      hist_bs2->Fill(Ds_btag_1,eventweight);
      hist_bt2->Fill(Dt_btag_1,eventweight);
      hist_bst_2d->Fill(Ds_btag_1,Dt_btag_1,eventweight);

      TopologicalVariables_new->Fill();

      count++;
      // if (Cut(ientry) < 0) continue;
   }

   cout << count << " events" << endl;
   outfile->Write();
   outfile->Close();

   cout << "delete file" << endl;
   delete outfile;
}



 
Double_t opt_disc::CalculateD_t(int Njet){

    double _Disc_tq=0.0;
    //  cout << acc_s << " " << acc_t << " " << c1 << " " << c2 << " " << c3 << " " << c4 << " " << c5 << " " << c6 << " " << c7 << endl; 

    if ( Njet == 2 ) { 
      if((P_tq + c1*P_wbb + c2*P_wcg + c3*P_wgg + c4*P_WW + c5*P_WZ + c6*P_ggg+ c7*P_lepjets_2jet)==0) {
	_Disc_tq=-1;
	cout << " Disc = - 1 "  << P_tq <<" "<< P_wbb<<" "<<P_wcg<<" "<<P_wgg<<" "<<P_WW<<" "<<P_WZ <<" "<<P_ggg<<" "<<P_lepjets_2jet<< endl;
      }
      
      else _Disc_tq = P_tq / ( P_tq + c1*P_wbb + c2*P_wcg + c3*P_wgg + c4*P_WW + c5*P_WZ + c6*P_ggg+ c7*P_lepjets_2jet);
      
    }
    else {
      if((P_tqb + c1*P_wbbg+ c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2+c4*P_wugg)==0) _Disc_tq=-1;
      else _Disc_tq = P_tqb / ( P_tqb + c1*P_wbbg + c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2+c4*P_wugg);
     }
    return _Disc_tq;
}


Double_t opt_disc::CalculateD_s(int Njet) {
  double _Disc_tb=0.;
  //  cout << acc_s << " " << acc_t << " " << c1 << " " << c2 << " " << c3 << " " << c4 << " " << c5 << " " << c6 << " " << c7 << endl; 

    if ( Njet == 2 ) { 
      if((P_tb + c1*P_wbb + c2*P_wcg + c3*P_wgg + c4*P_WW + c5*P_WZ + c6*P_ggg+ c7*P_lepjets_2jet)==0) _Disc_tb=-1;
      else _Disc_tb = P_tb / ( P_tb + c1*P_wbb + c2*P_wcg + c3*P_wgg + c4*P_WW + c5*P_WZ + c6*P_ggg+ c7*P_lepjets_2jet);
     
    }
    else {
      if((P_tbg + c1*P_wbbg+ c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2+c4*P_wugg)==0) {
	_Disc_tb=-1;
	cout << " Disc 3j = -1 "  <<" "<< P_tbg <<" "<< P_wbbg  <<" "<< P_lepjets_3jet_1 <<" "<< P_lepjets_3jet_2 <<endl; 
      }
      else _Disc_tb = P_tbg / ( P_tbg + c1*P_wbbg+ c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2+c4*P_wugg );
    }
    
    return _Disc_tb;
}


Double_t opt_disc::CalculateD_st(int Njet) {
    double _Disc_tb_tq=0.;
    //  cout << acc_s << " " << acc_t << " " << c1 << " " << c2 << " " << c3 << " " << c4 << " " << c5 << " " << c6 << " " << c7 << endl; 
    
    if ( Njet == 2 ) { 
      if((acc_s*P_tb + acc_t*P_tq + c1*P_wbb + c2*P_wcg + c3*P_wgg+ c4*P_WW + c5*P_WZ + c6*P_ggg+ c7*P_lepjets_2jet)==0) _Disc_tb_tq=-1;
      else _Disc_tb_tq = (acc_s*P_tb + acc_t*P_tq) / ( acc_s*P_tb + acc_t*P_tq + c1*P_wbb + c2*P_wcg + c3*P_wgg+ c4*P_WW + c5*P_WZ + c6*P_ggg+ c7*P_lepjets_2jet );
    }
    else {
      if( (acc_s*P_tbg + acc_t*P_tqb + c1*P_wbbg + c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2 +c4*P_wugg)==0) _Disc_tb_tq=-1;
      else _Disc_tb_tq = (acc_s*P_tbg + acc_t*P_tqb) / ( acc_s*P_tbg + acc_t*P_tqb + c1*P_wbbg + c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2+c4*P_wugg);
      }
    return _Disc_tb_tq;

}


Double_t opt_disc::CalculateD_btag_t(int Njet){

    double _Disc_tq=0.0;
    //  cout << acc_s << " " << acc_t << " " << c1 << " " << c2 << " " << c3 << " " << c4 << " " << c5 << " " << c6 << " " << c7 << endl; 

    if ( Njet == 2 ) { 
      if((P_tq + c1*P_wbb + c2*P_wcg + c3*P_wgg + c4*P_WW + c5*P_WZ + c6*P_ggg+ c7*P_lepjets_2jet)==0) _Disc_tq=-1;
      else _Disc_tq = bwgt*P_tq / ( bwgt*(P_tq + c1*P_wbb) + not_bwgt*(c2*P_wcg + c3*P_wgg + c4*P_WW + c5*P_WZ + c6*P_ggg) + c7*bwgt*P_lepjets_2jet);
      
    }
    else {
      if((P_tqb + c1*P_wbbg+ c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2+c4*P_wugg)==0) _Disc_tq=-1;
      else _Disc_tq = bwgt*P_tqb / ( bwgt*(P_tqb + c1*P_wbbg + c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2)+c4*not_bwgt*P_wugg);
     }
    return _Disc_tq;
}


Double_t opt_disc::CalculateD_btag_s(int Njet) {
  double _Disc_tb=0.;
  //  cout << acc_s << " " << acc_t << " " << c1 << " " << c2 << " " << c3 << " " << c4 << " " << c5 << " " << c6 << " " << c7 << endl; 

    if ( Njet == 2 ) { 
      if((P_tb + c1*P_wbb + c2*P_wcg + c3*P_wgg + c4*P_WW + c5*P_WZ + c6*P_ggg+ c7*P_lepjets_2jet)==0) _Disc_tb=-1;
      else _Disc_tb = bwgt*P_tb / ( bwgt*(P_tb + c1*P_wbb) + not_bwgt*(c2*P_wcg + c3*P_wgg + c4*P_WW + c5*P_WZ + c6*P_ggg)+ c7*bwgt*P_lepjets_2jet);
     
    }
    else {
      if((P_tbg + c1*P_wbbg+ c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2+c4*P_wugg)==0) _Disc_tb=-1;
      else _Disc_tb = bwgt*P_tbg / ( bwgt*(P_tbg + c1*P_wbbg+ c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2)+c4*not_bwgt*P_wugg );
    }
    
    return _Disc_tb;
}


Double_t opt_disc::CalculateD_btag_st(int Njet) {
    double _Disc_tb_tq=0.;
    //  cout << acc_s << " " << acc_t << " " << c1 << " " << c2 << " " << c3 << " " << c4 << " " << c5 << " " << c6 << " " << c7 << endl; 
    
    if ( Njet == 2 ) { 
      if((acc_s*P_tb + acc_t*P_tq + c1*P_wbb + c2*P_wcg + c3*P_wgg+ c4*P_WW + c5*P_WZ + c6*P_ggg+ c7*P_lepjets_2jet)==0) _Disc_tb_tq=-1;
      else _Disc_tb_tq = bwgt*(acc_s*P_tb + acc_t*P_tq) / ( bwgt*(acc_s*P_tb + acc_t*P_tq + c1*P_wbb) + not_bwgt*(c2*P_wcg + c3*P_wgg+ c4*P_WW + c5*P_WZ + c6*P_ggg)+ c7*not_bwgt*P_lepjets_2jet );
    }
    else {
      if( (acc_s*P_tbg + acc_t*P_tqb + c1*P_wbbg + c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2 +c4*P_wugg)==0) _Disc_tb_tq=-1;
      else _Disc_tb_tq = bwgt*(acc_s*P_tbg + acc_t*P_tqb) / ( bwgt*(acc_s*P_tbg + acc_t*P_tqb + c1*P_wbbg + c2*P_lepjets_3jet_1+c3*P_lepjets_3jet_2)+c4*not_bwgt*P_wugg);
      }
    return _Disc_tb_tq;

}