#include <iostream>
#include <fstream>
#include <sstream>
#include <stdio.h>
#include <map>
#include <string>
#include <set>
#include <iomanip>
#include <cmath>
#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <TLeaf.h>
#include <TLeafC.h>
#include <TH1F.h>
#include <TH2F.h>
#include <stdlib.h>

using namespace std;
int sizebegin=0;

map<string, int> _Check;
map<string, int> _CheckDatFile;
map<string, int> _MatchDatFile;
map<string, double> _TopoKey;
map<int, int> _CheckMatch;
bool Already;    



int main(int argc, char** argv) {

  Already = false;

  string inputfile;    
  int beginrange;
  string outputrootfilename;
  string rootfilename;
  int mc;
  int _NJets;
  string _Lepton;
  int _NTags;
  int btagging;
  string _Sig = "NoJES";
  if ( argc == 9 ) {
    stringstream j;
    j << argv[1];
    j >> _Lepton;
    cout << "Lepton: " << _Lepton << endl;
    stringstream k;
    k << argv[2];
    k >> _NJets;
    cout << "Jets: " << _NJets << endl;
    stringstream l;
    l << argv[3];
    l >> _NTags;
    cout << "Tags: " << _NTags << endl;
    stringstream m;
    m << argv[4];
    m >> rootfilename;
    cout << "TopovarTree File: " << rootfilename << endl;
    stringstream n;
    n << argv[5];
    n >> inputfile;
    cout << "Dat File: " << inputfile << endl; 
    stringstream q;
    q << argv[6];
    q >> outputrootfilename;
    cout << "Output RootFile: " << outputrootfilename << endl;
    stringstream r;
    r << argv[7];
    r >> mc;
    cout << "MC(0=data,qcd, 1=regular mc, 2=tb,tqb, 3=no mckey, 4=mckey but all the same, 5=match mc by lepton pt -> training sample 6=data,qcd with no mckey) " << mc << endl;
    stringstream s;
    s << argv[8];
    s >> btagging;
    cout << "Dat type(Inclusive=0, Btagging=1) " << btagging << endl;
    
  }
 
  else {
    cout << "You need to provide a number of jets, lepton, tags, rootfile, dat file, outputname" << endl;
    cout << argc << endl;
    for ( int i = 0; i < argc; i++ ) {
      cout << argv[i] << endl;
    }
    return(0);
  }


  double c1, c2, c3, c4, c5, c6, c7;
  double acc_s, acc_t;
    
  if ( _Lepton == "Electron" ) {
    if ( _NJets == 2 ) {
      if ( _NTags == 1 ) {
	c1 = 0.1;
	c2 = 0.4;
	c3 = 0.4;
	c4 = 0.1;
	acc_s = 7.0/18;
	acc_t = 11.0/18;
      }
      else if ( _NTags == 2 ) {
	c1 = 0.667;
	c2 = 0.0;
	c3 = 0.333;
	acc_s = 2.3/2.6;
	acc_t = 0.3/2.6;
      }
      else {
	cout << "Ntags not defined. Background fractions are not defined" << endl;
	exit(1);
      }
      }
    else if ( _NJets == 3 ) {
      if ( _NTags == 1 ) {
	c1 = 0.8;
	c2 = 0.2;
	acc_s = 3.0/9;
	acc_t = 6.0/9;
      }
      else if ( _NTags == 2 ) {
	c1 = 1.0;
	c2 = 0.0;
	acc_s = 1.1/1.9;
	acc_t = 0.8/1.9;
      }
      else {
	cout << "Ntags not defined. Background fractions are not defined" << endl;
	exit(1);
      }
      }
    else {
      cout << "_NJets not defined. Background fractions are not defined" << endl;
      exit(1);
    }
    }


  
  else if ( _Lepton == "Muon" ) {
    if ( _NJets == 2 ) {
      if ( _NTags == 1 ) {
	c1 = 0.04;
	c2 = 0.02;
	c3 = 0.88;
	c4 = .01;
	c5 = .01;
	c6 = .04;
	c7 = .001;
	acc_s = .36;
	acc_t = .64;
      }
      else if ( _NTags == 2 ) {
	c1 = 1.0;
	c2 = 0.0;
	c3 = 0.0;
	acc_s = 1.9/2.1;
	acc_t = 0.2/2.1;
      }
      else {
	cout << "Ntags not defined. Background fractions are not defined" << endl;
	exit(1);
      }
    }
    else if ( _NJets == 3 ) {
      if ( _NTags == 1 ) {
	c1 = 0.97;
	c2 = 0.015;
	c3 = 0.015;
	acc_s = .33;
	acc_t = .67;
      }
      else if ( _NTags == 2 ) {
	c1 = 1.0;
	c2 = 0.0;
	acc_s = 0.9/1.6;
	acc_t = 0.7/1.6;
      }
      else {
	cout << "Ntags not defined. Background fractions are not defined" << endl;
	exit(1);
      }
    }
    else {
      cout << "_NJets not defined. Background fractions are not defined" << endl;
      exit(1);
    }
  }
  else {
    cout << "Lepton not defined. Background fractions are not defined" << endl;
    exit(1);
  }
  
 
  ifstream input;
  input.open(inputfile.c_str());
  if ( ! input.is_open() ) {
    cout << "Could not open: " << inputfile << endl;
  }
  else {
    //    cout << "Opened: " << inputfile << endl;
  }
    

  TFile *RootFile = new TFile(rootfilename.c_str());
  if ( RootFile == NULL ) {
    cerr << "Could not open input file: " << rootfilename << endl;
    cerr << "Event.cpp" << endl;
  }	


  TTree *RootTree = (TTree*)RootFile->Get("TopologicalVariables");
  if ( RootTree == NULL ) {
    cerr << "Could not find root tree: TopologicalVariables" << endl;
  }
  
  int entries = RootTree->GetEntries();
  cout << "  Entries: " << entries << endl;
 
  //x cout << "Output Name: " << outputrootfilename << endl;


  TFile *file = new TFile(outputrootfilename.c_str(), "RECREATE");   
   
  TTree*  TopologicalVariables_new = RootTree->CloneTree(0);
  stringstream name1;
  name1 << "Discriminant";
  stringstream name2;
  name2 << "DiscriminantS";
  stringstream name3;
  name3 << "DiscriminantT";
  stringstream name4;
  name4 << "DiscriminantST";
  stringstream name5;
  name5 << "DiscriminantS_T";
  int bins = 9;
  double schannel[10]={-0.01, 0.22, 0.41, 0.57, 0.70, 0.79, 0.86, 0.91, 0.96, 1.01};
  double tchannel[10]={-0.01, 0.18, 0.39, 0.57, 0.70, 0.80, 0.88, 0.94, 0.98, 1.01};
  int st_bins = 18;
  double st_channel[19]={-0.01, 0.07, 0.16, 0.26, 0.36, 0.45, 0.54, 0.61, 0.67, 0.73, 0.78, 0.82, 0.86, 0.89, 0.92, 0.95, 0.97, 0.99, 1.01};  

  TH2F *hist = new TH2F(name1.str().c_str(), "Discriminant", bins, schannel, bins, tchannel);
  TH1F *hist_s = new TH1F(name2.str().c_str(), "Discriminant", bins, schannel);
  TH1F *hist_t = new TH1F(name3.str().c_str(), "Discriminant", bins, tchannel);
  TH1F *hist_st = new TH1F(name4.str().c_str(), "Discriminant", bins, tchannel);
  TH1F *hist_s_t = new TH1F(name5.str().c_str(), "Discriminant", st_bins, st_channel);

  double _Prob_tqg;
  double _Prob_tb;
  double _Prob_tbg;
  double _Prob_tq;
  double _Prob_tqb;
  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 _Disc_tb;
  double _Disc_tq;
  double _Disc_tb_tq;
  double _MissingEt;
  double _LeptonPt;
  double _Prob_wugg;
  double _jet2btagnn;
  double _jet1btagnn;
  double _jet3btagnn;
  


  TopologicalVariables_new->Branch("Disc_tb_btag_btag", &_Disc_tb, "_Disc_tb_btag/D");
  TopologicalVariables_new->Branch("Disc_tq_btag_btag", &_Disc_tq, "_Disc_tq_btag/D");
  TopologicalVariables_new->Branch("Disc_tb_tq_btag_btag", &_Disc_tb_tq, "_Disc_tb_tq_btag/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");
  
  int foundcount=0;
 
  int counter = -1;
  //for(int i=0; i < 50; i++) {
  for(int i=0; i < RootTree->GetEntries(); i++) {
    //cout << " hi" << endl;
    RootTree->GetEntry(i);
    //cout << RootTree->GetEntries() << endl;
    bool good = true;
      
    //consistency check...just in case
    //--------------------------------------------
    // Jet1
    //--------------------------------------------
    string jet1_pt_name = "Jet1Pt";
    double _jet1_pt = 0.0;
    TLeaf *jet1_pt_leaf = (TLeaf*)RootTree->GetLeaf(jet1_pt_name.c_str());
    if ( jet1_pt_leaf == NULL ) { cout << "Could not find " << jet1_pt_name << endl; good = false; }
    else { _jet1_pt = jet1_pt_leaf->GetValue(0); }

    // cout << _jet1_pt << " hello" << endl;
    
    //--------------------------------------------
    // Jet2
    //--------------------------------------------
    string jet2_pt_name = "Jet2Pt";
    double _jet2_pt = 0.0;
    TLeaf *jet2_pt_leaf = (TLeaf*)RootTree->GetLeaf(jet2_pt_name.c_str());
    if ( jet2_pt_leaf == NULL ) { cout << "Could not find " << jet2_pt_name << endl; good = false; }
    else { _jet2_pt = jet2_pt_leaf->GetValue(0); }

     
    //--------------------------------------------
    // Event Weights
    //--------------------------------------------
    string weight_total_name = "EventWeight";
    double _weight_total = -1.0;
    TLeaf *weight_total_leaf = (TLeaf*)RootTree->GetLeaf(weight_total_name.c_str());
    if ( weight_total_leaf == NULL ) { cout << "Could not find " << weight_total_name << endl; good = false; }
    else { _weight_total = weight_total_leaf->GetValue(0); }
      


    //--------------------------------------------
    // LeptonPt
    //--------------------------------------------
    string LeptonPt_name = "LeptonPt";
    TLeaf *LeptonPt_leaf = (TLeaf*)RootTree->GetLeaf(LeptonPt_name.c_str());
    if ( LeptonPt_leaf == NULL ) { cout << "Could not find " << LeptonPt_name << endl; good = false; }
    else { _LeptonPt = LeptonPt_leaf->GetValue(0); }


    //--------------------------------------------
    // Missing Et
    //--------------------------------------------
    string MissingEt_name = "METPt";
    TLeaf *MissingEt_leaf = (TLeaf*)RootTree->GetLeaf(MissingEt_name.c_str());
    if ( MissingEt_leaf == NULL ) { cout << "Could not find " << MissingEt_name << endl; good = false; }
    else { _MissingEt = MissingEt_leaf->GetValue(0); }

     //--------------------------------------------
    // Jet 1 BTag
    //--------------------------------------------
    string jet1btagnn_name = "Jet1BTagNN";
    TLeaf *jet1btagnn_leaf = (TLeaf*)RootTree->GetLeaf(jet1btagnn_name.c_str());
    if ( jet1btagnn_leaf == NULL ) { cout << "Could not find " << jet1btagnn_name << endl; good = false; }
    else { _jet1btagnn = jet1btagnn_leaf->GetValue(0); }

    //--------------------------------------------
    // Jet 2 BTag
    //--------------------------------------------
    string jet2btagnn_name = "Jet2BTagNN";
    TLeaf *jet2btagnn_leaf = (TLeaf*)RootTree->GetLeaf(jet2btagnn_name.c_str());
    if ( jet2btagnn_leaf == NULL ) { cout << "Could not find " << jet2btagnn_name << endl; good = false; }
    else { _jet2btagnn = jet2btagnn_leaf->GetValue(0); }

    //--------------------------------------------
    // Jet 3 BTag
    //--------------------------------------------
    string jet3btagnn_name = "Jet3BTagNN";
    TLeaf *jet3btagnn_leaf = (TLeaf*)RootTree->GetLeaf(jet3btagnn_name.c_str());
    if ( jet3btagnn_leaf == NULL ) { cout << "Could not find " << jet3btagnn_name << endl; good = false; }
    else { _jet3btagnn = jet3btagnn_leaf->GetValue(0); }

    int _jet1btag=0;
    int _jet2btag=0;
    int _jet3btag=0;

    //cout << _jet1btagnn << "-" <<  _jet2btagnn << "-" <<  _jet3btagnn << endl;

    if(_jet1btagnn > 0) _jet1btag=1;
    if(_jet2btagnn > 0) _jet2btag=1;
    if(_jet3btagnn > 0) _jet3btag=1;
    

    string _eventnumber_1 = "";
    if(mc==1 || mc ==2 || mc==3 || mc==4 || mc ==5) {
     //--------------------------------------------
    // Event Weights
    //--------------------------------------------
    string eventnumber_name = "MCkey";
    //char _eventnumber_1[];
    stringstream id;
    
    string _eventnumber_2 = "";
    //stringstream _eventnumber_1;
    TLeafC *eventnumber_leaf = dynamic_cast<TLeafC*>(RootTree->GetLeaf(eventnumber_name.c_str()));
    if ( eventnumber_leaf == NULL ) { cout << "Could not find " << eventnumber_name << endl; good = false; }
    else {_eventnumber_2 = eventnumber_leaf->GetValueString(); }
    //cout << "CHEKC : " << _eventnumber_2 << endl;
    
    if(mc == 3 || mc == 4) _eventnumber_2 = "nomckey" ;
    id << _eventnumber_2;
  
    if(mc==2) id << "-" << _jet1_pt << "-" << _jet2_pt;

    if(mc==5) id << "-" << int(_LeptonPt);    

    if(mc==3 || mc==4) {
      id << "-" << int(_jet1_pt) << "-" << int(_jet2_pt)<< "-" << int(_LeptonPt) << "-" <<int(_MissingEt);
    }
   
    if(btagging==1) {
      id << "-" <<_jet1btag <<_jet2btag;

      if(_NJets==3) id  <<_jet3btag ;
    }

    
    _eventnumber_1=id.str();
    
    }
    else {
    //--------------------------------------------
    // Jet1
    //--------------------------------------------
    string runnumber_name = "RunNumber";
    int runnumber = 0.0;
    TLeaf *runnumber_leaf = (TLeaf*)RootTree->GetLeaf(runnumber_name.c_str());
    if ( runnumber_leaf == NULL ) { cout << "Could not find " << runnumber_name << endl; good = false; }
    else { runnumber = runnumber_leaf->GetValue(0); }
    
    string eventname = "EventNumber";
    int eventnumber = 0.0;
    TLeaf *eventnumber_leaf = (TLeaf*)RootTree->GetLeaf(eventname.c_str());
    if ( eventnumber_leaf == NULL ) { cout << "Could not find " << eventname << endl; good = false; }
    else { eventnumber = eventnumber_leaf->GetValue(0); }
    
    stringstream id;
    id << eventnumber << "-" << runnumber;
    if(btagging==1) {
      id << "-" <<_jet1btag <<_jet2btag ;
      if(_NJets==3) id  <<_jet3btag ;
    }
     // <<"-" << int(_jet1_pt) << "-" << int(_jet2_pt) << "-" << int(_LeptonPt) << "-" << int(_MissingEt) ;
    _eventnumber_1=id.str();
    }
    
    //cout << " " <<i << " " << _eventnumber_1 << " event number " << endl;

    _Check[_eventnumber_1]=i;
    _CheckMatch[i]=0;


    /*
    //_TopoKey[_eventnumber_1]=_jet1_pt+_jet2_pt+_LeptonPt+_MissingEt;	
    //cout << _eventnumber_1 << endl;
    int topokey=(int)(1000*_jet1_pt)+ (int)(1000*_jet2_pt)+ (int)(1000*_LeptonPt)+ (int)(1000*_MissingEt)+ (int)(1000*_weight_total);
      
    int datkey =(int)(jet1pt*1000) + (int)(jet2pt*1000) + (int)(leptonpt*1000) +(int) (missinget*1000) + (int)(eventweight*1000);
    */

  }

  int xcount=0;
  cout << "_Check.size() is " << (int) _Check.size() << endl;

  /*
  for( map<string, int>::iterator ii=_Check.begin(); ii!=_Check.end(); ++ii)
    {
      xcount++;
      cout << xcount << " " << (*ii).first << ": " << (*ii).second << endl;
    }
  */

  int counter_mckey=0;
  int counter_nomckey=0;
  //while ( counter < 10) {
  int multiples=0;
  while ( input.eof() == false ) {
    counter++;

    string source;
    string lepton;
    string njets;
    string ntags;
    int run;
    int event;
    int eventno;
    int mcrandom1;
    int mcrandom2;
    double jet1pt;
    double jet2pt;
    double leptonpt;
    double missinget;
    double eventweight;
    double tb=-1.;
    double tbg=-1.;
    double tqperm1=-1.;
    double tqperm2=-1.;
    double tqb=-1.;
    double wbb=-1.;
    double wbbg=-1.;
    double wcgperm1=-1.;
    double wcgperm2=-1.;
    double WH105=-1.;
    double WH115=-1.;
    double WH125=-1.;
    double WH135=-1.;
    double WH145=-1.;
    double WH155=-1.;
    double WW=-1.;
    double WZ=-1.;
    double ggg=-1.;
    double wgg=-1.;
    double wugg= -1;
    double lepjets_2jet=-1.;
    double lepjets_3jet_1=-1.;
    double lepjets_3jet_2=-1.;
    double tqg=-1.;
    string dummy;
    string mckey="";
    string mckey1="";  
    input >> source >> lepton >> njets >> ntags >> run >> event >> eventno >> mcrandom1 >> mcrandom2 >> jet1pt >> jet2pt >> leptonpt >> missinget >> eventweight;


    if(mc==3 || mc ==4 ) mckey="nomckey";

    if(mc!=3 && mc!=6) input >> mckey;
     
    if(_NJets==2) {
       input >> tb >> tqperm1 >> tqperm2 >> wbb >> wcgperm1 >> wcgperm2 >> wgg >> WH105 >>  WH115 >>  WH125 >>  WH135 >>  WH145 >>  WH155 >> WW >> WZ >> lepjets_2jet >> ggg;
     }
    if(_NJets==3) {
       
       input >> tbg >> tqb >> tqg >> wbbg >> wugg >> lepjets_3jet_1 >> lepjets_3jet_2;
    }
    

    if ( source.size() == 0 ) { break; }
    
    if ( counter % 1000 == 0 && counter != 0 ) {
      //x      cout << "Read event: " << counter << "\t" << source << endl;
    }

    
    foundcount=0;
    int i = -1;
    if(mc==0 || mc==6) {
      stringstream id_1;
      id_1 << event << "-" << run ;
      if(btagging==1) id_1 << "-" << ntags; //<< "-" << int(jet1pt) << "-" << int(jet2pt) << "-" << int(leptonpt) << "-" << int(missinget);
      mckey1=id_1.str();
      //cout << mckey1 << " MCkey" << endl;
    }
    else {
      stringstream id_2;
      id_2 << mckey ;
      if (mc==2) id_2 << "-" << jet1pt << "-" << jet2pt;
      if (mc==3 || mc==4) id_2 << "-" << int(jet1pt) << "-" << int(jet2pt) << "-" << int(leptonpt) << "-" << int(missinget);
      if (mc==5) id_2 << "-" << int(leptonpt);
      if(btagging==1) id_2 <<"-"<<ntags; //<< "-" << int(leptonpt) << "-" << int(missinget);
      mckey1=id_2.str();
     

      //x      cout << " count key " << counter << " " << mckey1 << " " << ntags << endl;
    }
    _CheckDatFile[mckey1]=counter;
    map<string,int>::iterator iter = _Check.find(mckey1);
    if( iter != _Check.end() ) {
      foundcount++;	
      i=_Check.find(mckey1)->second;
      _CheckMatch[i]=1;
    }
    // In case there are multiple events in *.DAT then 
    // check to see if the event was already matched in written out.
    //if so, then mark it not found and move on
    map<string,int>::iterator itermatch = _MatchDatFile.find(mckey1);
    if( itermatch != _MatchDatFile.end() ) {
      foundcount=0;
      multiples++;
      //cout << " multiply matched evt in *.DAT " << _MatchDatFile.find(mckey1)->second << endl;
    }


    if(foundcount ==1) {
      _MatchDatFile[mckey1]=i;
      RootTree->GetEntry(i);
      _Prob_tb = tb;
      _Prob_tbg = tbg;
      _Prob_tq =  tqperm1+tqperm2;
      _Prob_tqb = tqb;
      _Prob_tqg = tqg;
      _Prob_wbb = wbb;
      _Prob_wcg = wcgperm1+wcgperm2;
      _Prob_wgg = wgg;
      _Prob_wbbg = wbbg;
      _Prob_WW = WW;
      _Prob_WZ= WZ;
      _Prob_ggg= ggg;
      _Prob_lepjets_2jet = lepjets_2jet;
      _Prob_lepjets_3jet_1 = lepjets_3jet_1;
      _Prob_lepjets_3jet_2 = lepjets_3jet_2;
      _Prob_wugg = wugg;
      if(_Prob_tb < 0) _Prob_tb =0.0;
      if(_Prob_tbg < 0) _Prob_tbg =0.0;
      if(_Prob_tqg < 0) _Prob_tqg =0.0;
      if(_Prob_tq < 0) _Prob_tq =0.0;
      if(_Prob_tqb < 0) _Prob_tqb =0.0;
      if(_Prob_wbb < 0) _Prob_wbb =0.0;
      if(_Prob_wcg < 0) _Prob_wcg =0.0;
      if(_Prob_wgg < 0) _Prob_wgg =0.0;
      if(_Prob_wbbg < 0) _Prob_wbbg =0.0;
      if(_Prob_WW < 0) _Prob_WW =0.0;
      if(_Prob_WZ < 0) _Prob_WZ =0.0;
      if(_Prob_ggg < 0) _Prob_ggg =0.0;
      if(_Prob_lepjets_2jet < 0) _Prob_lepjets_2jet =0.0;
      if(_Prob_lepjets_3jet_1 < 0) _Prob_lepjets_3jet_1 =0.0;
      if(_Prob_lepjets_3jet_2 < 0) _Prob_lepjets_3jet_2 =0.0;
      if( _Prob_wugg <0) _Prob_wugg =0.0;
      if ( _NJets == 2 ) {
	_Disc_tb = _Prob_tb / ( _Prob_tb + c1*_Prob_wbb + c2*_Prob_wcg + c3*_Prob_wgg + c4*_Prob_WW + c5*_Prob_WZ + c6*_Prob_ggg+ c7*_Prob_lepjets_2jet);
	_Disc_tq = _Prob_tq / ( _Prob_tq + c1*_Prob_wbb + c2*_Prob_wcg + c3*_Prob_wgg + c4*_Prob_WW + c5*_Prob_WZ + c6*_Prob_ggg+ c7*_Prob_lepjets_2jet);
	_Disc_tb_tq = (acc_s*_Prob_tb + acc_t*_Prob_tq) / ( acc_s*_Prob_tb + acc_t*_Prob_tq + c1*_Prob_wbb + c2*_Prob_wcg + c3*_Prob_wgg+ c4*_Prob_WW + c5*_Prob_WZ + c6*_Prob_ggg+ c7*_Prob_lepjets_2jet );
      }
      else {
	_Disc_tb = _Prob_tbg / ( _Prob_tbg + c1*_Prob_wbbg+ c2*_Prob_lepjets_3jet_1+c3*_Prob_lepjets_3jet_2 +c4*_Prob_wugg);
	_Disc_tq = (.5*_Prob_tqb+.5*_Prob_tqg) / ( .5*_Prob_tqb+.5*_Prob_tqg + c1*_Prob_wbbg + c2*_Prob_lepjets_3jet_1+c3*_Prob_lepjets_3jet_2+c4*_Prob_wugg);
	_Disc_tb_tq = (acc_s*_Prob_tbg + acc_t*(.5*_Prob_tqb+.5*_Prob_tqg)) / ( acc_s*_Prob_tbg + acc_t*( .5*_Prob_tqb+.5*_Prob_tqg)+ c1*_Prob_wbbg + c2*_Prob_lepjets_3jet_1+c3*_Prob_lepjets_3jet_2+c4*_Prob_wugg);
      }

      if ( _Prob_tb == 0.0 ||
	   _Prob_tq == 0.0 ||
	   _Prob_wbb == 0.0 ||
	   _Prob_wcg == 0.0 ||
	   _Prob_wgg == 0.05 ) {
	//	  cout << "Something bad happened" << endl;
	//	  cout << id << endl;
	//	  exit(1);
      }
      if ( _Disc_tb != _Disc_tb ) { _Disc_tb = -1.0; }
      if ( _Disc_tq != _Disc_tq ) { _Disc_tq = -1.0; }
      if ( _Disc_tb_tq != _Disc_tb_tq ) { _Disc_tb_tq = -1.0; }
      //tree->Fill();
      hist->Fill(_Disc_tb, _Disc_tq, eventweight);
      hist_s->Fill(_Disc_tb, eventweight);
      hist_t->Fill(_Disc_tq, eventweight);
      hist_st->Fill(_Disc_tb_tq, eventweight);
      hist_s_t->Fill(_Disc_tb_tq, eventweight);
      TopologicalVariables_new->Fill();
      //	cout << _Disc_tb << "\t" << _Disc_tq << "\t" << hist->GetEntries() << endl;
      }
    else {	
     
      //cout << "Problem: no probability " << endl;
      //	cout << outputrootfilename.str() << "\t" << _NoProbWeight <<
      //	"\t" << id << endl;
      //	exit(1);
    }
    //RootTree->GetEntry(entry);
    
  }
  
  for(int i=0; i< _CheckMatch.size() ; i++) {
    //map<string,int>::iterator itercheckmatch = _CheckMatch.find(i);
    if( _CheckMatch.find(i)->second == 0 ) cout << i << " entry in root file needs to be looked into" << endl;
  }

  //if( TopologicalVariables_new->GetEntries()!= counter) cout << "!!!!!!!!!Problem!!!!!!!!!!!!! Not all matched" <<endl;
  //cout << " No mckey .dat: " << counter_nomckey << " MCkey in .dat: " << counter_mckey << endl;
  cout << " numEVTS in DatFile " << counter << " _CheckDatFile.size() is " << (int) _CheckDatFile.size() << " and matched " << _MatchDatFile.size() << " and multiples " << multiples << endl ;   
  cout << " numEVTS in input RootFile " << RootTree->GetEntries() << " num w/o duplicate mckey " << _Check.size() << endl;
  cout << " numEVTS in OUTPUT RootFile " << TopologicalVariables_new->GetEntries() << endl;
  
  if( TopologicalVariables_new->GetEntries()!= _Check.size()) cout << inputfile << " !!!!!!!!!Problem!!!!!!!!!!!!! Numbers DO NOT  match " << endl;


  file->Write();
  file->Close();
 
  cout << "delete file" << endl;
  delete file;
  delete RootFile;

  cout << "done" << endl;

  return( 0 );
}