#include "AnglesUtil.hpp"

#include "TClonesArray.h"
#include <TRandom3.h>
void IntHist(TH1* hist, TH1* hist_int)
{
  for (int bin=0; bin<hist->GetNbinsX(); bin++)
    {
      //cout<<hist->Integral(0,hist->GetNbinsX())<<"  "<<hist->GetNbinsX()<<endl;;
      hist_int->SetBinContent(bin,hist->Integral(bin,hist->GetNbinsX())/hist->Integral(0,hist->GetNbinsX()));
    }
}


void SetI_SM_mlpTop_type3_tauQCD(Int_t ntrain=100) {


  // For a LEP search for invisible Higgs boson, a neural network 
  // was used to separate the signal from the background passing 
  // some selection cuts. Here is a simplified version of this network, 
  // taking into account only WW events.
  //    Author: Christophe Delaere

  gSystem->CompileMacro("TopTree_classes/TopTauTreeClasses.C","k");
  
  
  if (!gROOT->GetClass("TMultiLayerPerceptron")) {
    //gSystem->Load("./libHist.so.5.15");
    gSystem->Load("libMLP.so");      
  }
  
  TClonesArray*  taus  = new TClonesArray("TopTau",1);  
  TClonesArray*  jets  = new TClonesArray("TopJet",1);
  
  

  // Prepare inputs
  // The 2 trees are merged into one, and a "type" branch,
  // equal to 1 for the signal and 0 for the background is added.

  // I N P U T S
  TChain *background = new TChain("TopTauTree");

   background->Add("/work/taichi-clued0/dmenezes/DATA_VCtaujets/vcdata*.root");
   TFile inputsig("/work/mengo-clued0/dmenezes/output_MC/NOMETSIG/TTau_ttbarFullp20M170AlpSMTauJetsMC.root");
  
  // I N I T I A L I Z E  C U T  V A R I A B L E S
  Float_t TAUET      = 10.0;
  Float_t TAUETMAX   = 150.0;
  Float_t TAUETA     = 2.50;
  Float_t TAUNN      = 0.95;
  Float_t TAUNNMIN   = 0.90;
  Float_t TAUNNelec  = -1.00;
  Float_t MINCRACK   = -100.0;
  Float_t MAXCRACK   = 100.0;    
  Float_t JETPT      = 25.0;
  Float_t JETPTMAX   = 200.0;
  Float_t JETETA     = 2.50;
  Float_t METL     = 4.0;
  Int_t   NJETS      = -3;
  Int_t   NTAG       = 1;

  TTree *signal      = (TTree *) inputsig.Get("TopTauTree");
  TFile *output      = new TFile("Hplus_OUTPUTS/TypeIII/TrainingOutPut/SetI_SM_type3Set15-40.root","RECREATE","select",5);
  output->cd();
  TTree *simu        = new TTree("TrainTree", "Tree for training NN");
  

  // I N I T I A L I Z E  B R A N C H  V A R I A B L E S
  Int_t type,evtno;
  Float_t nn, taubestNN, weight;
  Float_t aplan, btag_mcweight, cent, costhetastar, dphiemet, dphitaumet,dphitaumet_mc,dphitaumet_QCD, electronEta, electronPhi, electronPt,geommeanPt, ht, ktmin, ktminp, lcorr_mcweight, met, metl, metx, mety, min2jetmass, mtemet, mttaumet, mttaumet_mc, mttaumet_QCD, spher, sqrts, topmass, topmassl, trig_mcweight, trig_mcweight_minus, trig_mcweight_plus, weightedetarms, wmass, wmassl, McMet, A, sigma, l1, l2, met_noTauJes, metl_noTauJes, met_phi, dphitaumet_noTauJes, met_perp_noTauJes, met_colin_noTauJes;
  Double_t lum_mcweight, lum_reweight, PVz_reweight;
  Double_t trigeff_0btag, trigeff_1btag, trigeff_2btag, trigeff_3btag;
  Double_t tau_nnout_corr_p20, tau_track_corr_p20,bFragWeight;
  
  // S i g n a l  t t b a r   - >  T a u
  signal->SetBranchAddress("TopTau",  &taus);
  signal->SetBranchAddress("TopJet",  &jets); 
  signal->SetBranchAddress("aplan",  &aplan);
  signal->SetBranchAddress("btag_mcweight", &btag_mcweight);
  signal->SetBranchAddress("lum_reweight", &lum_reweight);
  signal->SetBranchAddress("PVz_reweight", &PVz_reweight);
  signal->SetBranchAddress("cent",  &cent);
  signal->SetBranchAddress("costhetastar", &costhetastar);
  signal->SetBranchAddress("dphiemet",  &dphiemet);
  signal->SetBranchAddress("dphitaumet",  &dphitaumet);
  signal->SetBranchAddress("electronEta",  &electronEta);
  signal->SetBranchAddress("electronPhi",  &electronPhi);
  signal->SetBranchAddress("electronPt",  &electronPt);
  signal->SetBranchAddress("geommeanPt",  &geommeanPt);
  signal->SetBranchAddress("ht",  &ht);
  signal->SetBranchAddress("ktmin",  &ktmin);
  signal->SetBranchAddress("ktminp",  &ktminp);
  signal->SetBranchAddress("lcorr_mcweight", &lcorr_mcweight);
  signal->SetBranchAddress("lum_mcweight", &lum_mcweight);
  signal->SetBranchAddress("met",  &met);
  signal->SetBranchAddress("metl",  &metl);
  signal->SetBranchAddress("metx",  &metx);
  signal->SetBranchAddress("mety",  &mety);
  signal->SetBranchAddress("min2jetmass",  &min2jetmass);
  signal->SetBranchAddress("mtemet",  &mtemet);
  signal->SetBranchAddress("mttaumet",  &mttaumet);
  signal->SetBranchAddress("spher",  &spher);
  signal->SetBranchAddress("sqrts",  &sqrts);
  signal->SetBranchAddress("topmass",  &topmass);
  signal->SetBranchAddress("topmassl",  &topmassl);
  signal->SetBranchAddress("trig_mcweight", &trig_mcweight);
  signal->SetBranchAddress("trig_mcweight_minus", &trig_mcweight_plus);
  signal->SetBranchAddress("trig_mcweight_plus", &trig_mcweight_plus);
  signal->SetBranchAddress("weightedetarms", &weightedetarms);
  signal->SetBranchAddress("wmass", &wmass);
  signal->SetBranchAddress("wmassl", &wmassl);
  signal->SetBranchAddress("A", &A);
  signal->SetBranchAddress("l1", &l1);
  signal->SetBranchAddress("l2", &l2);	
  signal->SetBranchAddress("sigma", &sigma);
  signal->SetBranchAddress("trigeff_0btag", &trigeff_0btag);
  signal->SetBranchAddress("trigeff_1btag", &trigeff_1btag);
  signal->SetBranchAddress("trigeff_2btag", &trigeff_2btag);
  signal->SetBranchAddress("trigeff_3btag", &trigeff_3btag);
  signal->SetBranchAddress("tau_nnout_corr_p20", &tau_nnout_corr_p20);
  signal->SetBranchAddress("tau_track_corr_p20", &tau_track_corr_p20);
  signal->SetBranchAddress("bFragWeight", &bFragWeight);
  signal->SetBranchAddress("evtno", &evtno);
  // 4  J e t  S k i m  D A T A (B a c k g r o u n d)
  background->SetBranchAddress("TopTau",  &taus);
  background->SetBranchAddress("TopJet",  &jets); 
  background->SetBranchAddress("aplan",  &aplan);
  background->SetBranchAddress("btag_mcweight", &btag_mcweight);
  background->SetBranchAddress("lum_reweight", &lum_reweight);
  background->SetBranchAddress("PVz_reweight", &PVz_reweight);
  background->SetBranchAddress("cent",  &cent);
  background->SetBranchAddress("costhetastar", &costhetastar);
  background->SetBranchAddress("dphiemet",  &dphiemet);
  background->SetBranchAddress("dphitaumet",  &dphitaumet);
  background->SetBranchAddress("electronEta",  &electronEta);
  background->SetBranchAddress("electronPhi",  &electronPhi);
  background->SetBranchAddress("electronPt",  &electronPt);
  background->SetBranchAddress("geommeanPt",  &geommeanPt);
  background->SetBranchAddress("ht",  &ht);
  background->SetBranchAddress("ktmin",  &ktmin);
  background->SetBranchAddress("ktminp",  &ktminp);
  background->SetBranchAddress("lcorr_mcweight", &lcorr_mcweight);
  background->SetBranchAddress("lum_mcweight", &lum_mcweight);
  background->SetBranchAddress("met",  &met);
  background->SetBranchAddress("metl",  &metl);
  background->SetBranchAddress("metx",  &metx);
  background->SetBranchAddress("mety",  &mety);
  background->SetBranchAddress("min2jetmass",  &min2jetmass);
  background->SetBranchAddress("mtemet",  &mtemet);
  background->SetBranchAddress("mttaumet",  &mttaumet);
  background->SetBranchAddress("spher",  &spher);
  background->SetBranchAddress("sqrts",  &sqrts);
  background->SetBranchAddress("topmass",  &topmass);
  background->SetBranchAddress("topmassl",  &topmassl);
  background->SetBranchAddress("trig_mcweight", &trig_mcweight);
  background->SetBranchAddress("trig_mcweight_minus", &trig_mcweight_plus);
  background->SetBranchAddress("trig_mcweight_plus", &trig_mcweight_plus);
  background->SetBranchAddress("weightedetarms", &weightedetarms);
  background->SetBranchAddress("wmass", &wmass);
  background->SetBranchAddress("wmassl", &wmassl);
  background->SetBranchAddress("A", &A);
  background->SetBranchAddress("l1", &l1);
  background->SetBranchAddress("l2", &l2);	
  background->SetBranchAddress("sigma", &sigma);
//  background->SetBranchAddress("trigeff_0btag", &trigeff_0btag);
//  background->SetBranchAddress("trigeff_1btag", &trigeff_1btag);
//  background->SetBranchAddress("trigeff_2btag", &trigeff_2btag);
//  background->SetBranchAddress("trigeff_3btag", &trigeff_3btag);
  // S e t t i n g  U p  S i m u l a t i o n  B r a n c h e s   
  simu->Branch("metl",         &metl,         "metl/F");
  simu->Branch("mttaumet",     &mttaumet,     "mttaumet/F");
  simu->Branch("topmassl",     &topmassl,     "topmassl/F");
  simu->Branch("spher",        &spher,        "spher/F");
  simu->Branch("aplan",        &aplan,        "aplan/F");
  simu->Branch("cent",         &cent,         "cent/F");
  simu->Branch("sqrts",        &sqrts,        "sqrts/F");
  simu->Branch("ht",           &ht,           "ht/F");
  simu->Branch("ktminp",       &ktminp,       "ktminp/F");
  simu->Branch("costhetastar", &costhetastar, "costhetastar/F");
  simu->Branch("type",         &type,         "type/I");
  // simu->Branch("taubestNN",    &taubestNN,    "taubestNN/F");
  simu->Branch("metl_noTauJes",    &metl_noTauJes,    "metl_noTauJes/F");
  simu->Branch("weight",       &weight,       "weight/F");
  simu->Branch("btag_mcweight",       &btag_mcweight,       "btag_mcweight/F");
  simu->Branch("lum_mcweight",       &lum_mcweight,       "lum_mcweight/D");
  simu->Branch("trig_mcweight",       &trig_mcweight,       "trig_mcweight/F");
  simu->Branch("btag_mcweight",       &bFragWeight,       "bFragWeight/D");
  simu->Branch("lum_mcweight",       &PVz_reweight,       "PVz_reweight/D");
  simu->Branch("trig_mcweight",       &lum_reweight,       "lum_reweight/D");
  simu->Branch("evtno",       &evtno,       "evtno/I");
  
  type            = 1; // meaning signal 
  Float_t num_sig = 0.0;
  Float_t num_sig_true  = 0.0;
  Float_t num_bg  = 0.0;
  Int_t i;
  
  cout << "signal_tau+jets->GetEntries() = " << signal->GetEntries() << endl;

//  for(i = 0; i < signal->GetEntries()-2700; i++){
  for(i = 0; i < signal->GetEntries(); i++){
    signal->GetEntry(i);
    taus->Sort();
    
/////////////////////////////////////SIGNAL SELECTION  STARTS HERE///////////////////////////////////////////////////////////////
    
//    if (met>200 || metl < METL) continue;
    if (met>200) continue;
    Int_t sel_taus = 0;
    Int_t sel_taus_QCD = 0;
    Int_t sel_taus_typeOP = 0;
    Int_t sel_jets = 0;
    
    TopTau* tau0;
    TopTau* tau0QCD;
    
    Float_t sum_taus_QCD = 0;
    
    Int_t num_tau = taus->GetLast()+1;
    if(num_tau>0){

      for(Int_t itau = 0; itau < num_tau; itau++){
	TopTau* tau  = (TopTau*)taus->At(itau);
	// cout << " Tau "<<itau<<"; NN = "<<tau->_nnout<<"; Type = "<<tau->_tautype<<endl;
	// S e l e c t i n g  S i g n a l  T a u s
	if(tau->_Et>TAUET && tau->_Et<TAUETMAX && tau->_nnout>TAUNN && fabs(tau->_etad)<TAUETA && (tau->_tautype == 3)  && tau->_phicrack > MINCRACK && tau->_phicrack < MAXCRACK ){
	  if(sel_taus==0) tau0 = tau;
	  if(tau0->_tautype == 3){
	    // cout<<  " --------------- Sel Tau ------------ "<<endl;
	    // cout << " Tau "<<itau<<"; NN = "<<tau->_nnout<<"; Type = "<<tau->_tautype<<endl;
	    // cout<<  " --------------- Sel Tau ------------ "<<endl;
	    sel_taus++;
	  }
	}
	
	// S e l e c t i n g  Q C D  T a u s
	if(tau->_Et>TAUET && tau->_Et<TAUETMAX && tau->_nnout>0.3 && tau->_nnout<0.7 && fabs(tau->_eta)<TAUETA && (tau->_tautype == 3)  && tau->_phicrack > MINCRACK && tau->_phicrack < MAXCRACK ){
	  sum_taus_QCD+=tau->_nnout;
	  sel_taus_QCD++;
	}          
      }
      TRandom3 *RandomNumber = new TRandom3();

      Float_t R = RandomNumber->Rndm()*sum_taus_QCD;
      sum_taus_QCD = 0;   
      for(Int_t itau = 0; itau < num_tau; itau++){
	TopTau* tau  = (TopTau*)taus->At(itau);
	if(tau->_Et>TAUET && tau->_Et<TAUETMAX && tau->_nnout>0.3 && tau->_nnout<0.70 && fabs(tau->_eta)<TAUETA && (tau->_tautype == 3)  && tau->_phicrack > MINCRACK && tau->_phicrack < MAXCRACK ){
	  sum_taus_QCD+=tau->_nnout;
	  if(sum_taus_QCD>R) {tau0QCD = tau; break;}
	}         
	
      }
    }
    delete  RandomNumber; 
    
    // J e t  S e l e c t i o n  S t a r t s  h e r e
    Float_t btag_mcweight_0exec = 0;
    Float_t btag_mcweight_1exec = 0;
    Float_t btag_mcweight_2 = 0;
    Float_t btag_mcweight_2exec = 0;
    Float_t btag_mcweight_3 = 0;
    
    Float_t btag_mcweight = 1;  
    Int_t num_jet   = jets->GetLast()+1;
    for(Int_t ijet = 0; ijet < num_jet; ijet++){
      TopJet* jet   = (TopJet*)jets->At(ijet);
      Float_t pt    = sqrt(jet->_px*jet->_px+jet->_py*jet->_py);
      
      if(pt>JETPT  && pt<JETPTMAX && fabs(jet->_etad/10.0)<JETETA && jet->_dR>0.5){
	 btag_mcweight *= 1 - jet->_data_trf_nn;
	 sel_jets++;

   Float_t btag_mcweight_ijet_exec = jet->_data_trf_nn;      
      
    for (Int_t jjet = 0; jjet < num_jet; jjet++){
      TopJet* jet2 = (TopJet*)jets->At(jjet);
       Float_t pt = sqrt(jet2->_px*jet2->_px+jet2->_py*jet2->_py);
	
       if(pt>JETPT  && pt<JETPTMAX && fabs(jet2->_etad/10.0)<JETETA && jet2->_dR>0.5 && jjet!=ijet){
        
	        btag_mcweight_ijet_exec *= 1 - jet2->_data_trf_nn;	

       		Float_t btag_mcweight_ijet_or_jjet_exec = jet->_data_trf_nn*jet2->_data_trf_nn;
	  
                
	   	for (Int_t tjet = 0; tjet < num_jet; tjet++){
             		TopJet* jet3 = (TopJet*)jets->At(tjet);
               		Float_t pt = sqrt(jet3->_px*jet3->_px+jet3->_py*jet3->_py);
			                          
			if(pt>JETPT  && pt<JETPTMAX && fabs(jet3->_etad/10.0)<JETETA && jet3->_dR>0.5 && tjet!=ijet && tjet!=jjet)
			{
				btag_mcweight_ijet_or_jjet_exec *= 1 - jet3->_data_trf_nn;
				//cout<<"ijet = "<<ijet<<" jjet = "<<jjet<<" tjet = "<<tjet<<endl;
				//cout<<"jet->_data_trf_nn = "<<jet->_data_trf_nn<<" jet2->_data_trf_nn = "<<jet2->_data_trf_nn<<
				//" jet3->_data_trf_nn = "<<jet3->_data_trf_nn<<endl;
				//cout<<"btag_mcweight_ijet_or_jjet_execA = "<<btag_mcweight_ijet_or_jjet_exec<<endl;
			}	
		                                 
		} 
	 	//cout<<"btag_mcweight_ijet_or_jjet_execB = "<<btag_mcweight_ijet_or_jjet_exec<<endl;
         	btag_mcweight_2exec += btag_mcweight_ijet_or_jjet_exec;
	}
					    
     }
   
     btag_mcweight_1exec += btag_mcweight_ijet_exec;

      }
    }

///    btag_mcweight = 1 - btag_mcweight; the last one I commented

    btag_mcweight_2exec = btag_mcweight_2exec / 2;
    btag_mcweight = 1 - btag_mcweight;
    btag_mcweight_0exec = 1 - btag_mcweight;
    btag_mcweight_2 = btag_mcweight - btag_mcweight_1exec;
    btag_mcweight_3 = btag_mcweight_2 - btag_mcweight_2exec;

	trig_mcweight = trigeff_0btag*btag_mcweight_0exec + trigeff_1btag*btag_mcweight_1exec +  trigeff_2btag*btag_mcweight_2exec + trigeff_3btag*btag_mcweight_3;  


    //*******************************************************************

    // T h i s  i s  f o r  p y t h i a  s a m p l e s  o n l y
//    lum_mcweight = 1;
    if (sel_taus<=0 && sel_taus_QCD<=0) continue;
    if (sel_taus<=0) tau0 = tau0QCD;
    
    // F i n d   T h e  D e l t a  R  b e tw e e n  T a u s (tau cone of 0.3)  a n d  J e t s (jet cone of 0.5)
    Float_t dR = 100000;
    for (Int_t ijet = 0; ijet < num_jet; ijet++){
      TopJet* jet   = (TopJet*)jets->At(ijet);
      Float_t pt    = sqrt(jet->_px*jet->_px+jet->_py*jet->_py);
      
      if(pt/jet->_jes>JETPT  && pt/jet->_jes<JETPTMAX && fabs(jet->_etad/10.0)<JETETA)
	{
	  Float_t dRtmp = kinem::delta_R(tau0->_etad,tau0->_phi,jet->_etad/10,jet->_phi);
	  if(dRtmp<dR) {dR = dRtmp; taujet=jet;}
	} 
    } 


    Int_t Tau_Type = tau0->_tautype; 
    Float_t tauEta = tau0->_etad;
    Float_t tauPhi = tau0->_phi;
    Float_t tauPt  = tau0->_Et;
    met_phi = kinem::phi(metx,mety); 
    dphitaumet    = kinem::delta_phi(tau0->_phi, met_phi);
    mttaumet      = TMath::Sqrt(2*tau0->_Et*met*(1-TMath::Cos(dphitaumet)));
    TopTau* tau00 = (TopTau*)taus->At(0);
    if(met>0) ktminp = ktminp*(tau00->_Et+met)/met; else ktminp = -1;
    
    // F i l l i n g  H i s t o g r a m s F o r  T i g h t  T a u s
    
    if(sel_taus>0 && sel_jets>NJETS && metl>METL && (evtno%3==0)){
//      weight=lum_mcweight*trig_mcweight*btag_mcweight*lum_reweight*PVz_reweight*bFragWeight*tau_nnout_corr_p20*tau_track_corr_p20*141;     
//      num_sig+=lum_mcweight*trig_mcweight*btag_mcweight*lum_reweight*PVz_reweight*bFragWeight*tau_nnout_corr_p20*tau_track_corr_p20*141;
      weight=lum_mcweight*trig_mcweight*btag_mcweight*lum_reweight*PVz_reweight*bFragWeight*140;     
      num_sig+=lum_mcweight*trig_mcweight*btag_mcweight*lum_reweight*PVz_reweight*bFragWeight*140;
      num_sig_true++;
      simu->Fill();


    }
  }

/////////////////////////////////////SIGNAL SELECTION ENDS HERE///////////////////////////////////////////////////////////////

  
  type = 0;
  for(i = 0; i <= 600000; i++){
    background->GetEntry(i);
    taus->Sort();
    if(i%10000==0) cout<<i<<" / "<<background->GetEntries()<<endl;
//    if (met>200 || metl < METL) continue;


/////////////////////////////////////BACKGROUND SELECTION  STARTS HERE///////////////////////////////////////////////////////////////

    if (met>200) continue;
    
    Int_t sel_taus = 0;
    Int_t sel_taus_QCD = 0;
    Int_t sel_jets = 0;
    Int_t sel_jets_multijet = 0;
    Int_t sel_jetstag = 0;
    Int_t sel_jetstaggable = 0;
    Int_t sel_taus_typeOP = 0;
    TopJet* taujet;
    
    TopTau* tau0;
    TopTau* tau0QCD;
    
    Float_t sum_taus_QCD = 0;
    
    Int_t num_tau = taus->GetLast()+1;
    if(num_tau>0){
      for(Int_t itau = 0; itau < num_tau; itau++){
	TopTau* tau  = (TopTau*)taus->At(itau);
	// cout << " Tau "<<itau<<"; NN = "<<tau->_nnout<<"; Type = "<<tau->_tautype<<endl;
	// S e l e c t i n g  S i g n a l  T a u s
	if(tau->_Et>TAUET && tau->_Et<TAUETMAX && tau->_nnout>TAUNN && fabs(tau->_etad)<TAUETA && (tau->_tautype == 3)  && tau->_phicrack > MINCRACK && tau->_phicrack < MAXCRACK ){
	  if(sel_taus==0) tau0 = tau;
	  if(tau0->_tautype == 3){
	    // cout<<  " --------------- Sel Tau ------------ "<<endl;
	    // cout << " Tau "<<itau<<"; NN = "<<tau->_nnout<<"; Type = "<<tau->_tautype<<endl;
	    // cout<<  " --------------- Sel Tau ------------ "<<endl;
	    sel_taus++;
	  }
	}
	
	// S e l e c t i n g  Q C D  T a u s
	if(tau->_Et>TAUET && tau->_Et<TAUETMAX && tau->_nnout>0.3 && tau->_nnout<0.70 && fabs(tau->_eta)<TAUETA && (tau->_tautype == 3)  && tau->_phicrack > MINCRACK && tau->_phicrack < MAXCRACK ){
	  sum_taus_QCD+=tau->_nnout;
	  sel_taus_QCD++;
	}          
      }
      TRandom3 *RandomNumber = new TRandom3();

      Float_t R = RandomNumber->Rndm()*sum_taus_QCD;
      sum_taus_QCD = 0;   
      for(Int_t itau = 0; itau < num_tau; itau++){
	TopTau* tau  = (TopTau*)taus->At(itau);
	if(tau->_Et>TAUET && tau->_Et<TAUETMAX && tau->_nnout>0.3 && tau->_nnout<0.70 && fabs(tau->_eta)<TAUETA && (tau->_tautype == 3)  && tau->_phicrack > MINCRACK && tau->_phicrack < MAXCRACK ){
	  sum_taus_QCD+=tau->_nnout;
	  if(sum_taus_QCD>R) {tau0QCD = tau; break;}
	}         
	
      }
    }
    delete  RandomNumber;  
    // J e t  s e l e c t i o n  s t a r t s  h e r e 
    Int_t sel_jetstag = 0;    
    Int_t num_jet     = jets->GetLast()+1;
    for (Int_t ijet = 0; ijet < num_jet; ijet++){
      TopJet* jet = (TopJet*)jets->At(ijet);
      Float_t pt  = sqrt(jet->_px*jet->_px+jet->_py*jet->_py);
       
      if(pt>JETPT  && pt<JETPTMAX && fabs(jet->_etad/10.0)<JETETA && jet->_dR>0.5 && jet->_tagged_nn){
	sel_jetstag++;
      }
    }
    
    if (sel_taus<=0 && sel_taus_QCD<=0) continue;
    if (sel_taus<=0) tau0 = tau0QCD; // if there is no "good" taus - pick whichever we have, meaning the QCD tau!
    // F i n d   T h e  D e l t a  R  b e tw e e n  T a u s (tau cone of 0.3)  a n d  J e t s (jet cone of 0.5)
    Float_t dR = 100000;
    for (Int_t ijet = 0; ijet < num_jet; ijet++){
      TopJet* jet   = (TopJet*)jets->At(ijet);
      Float_t pt    = sqrt(jet->_px*jet->_px+jet->_py*jet->_py);
      
      if(pt/jet->_jes>JETPT  && pt/jet->_jes<JETPTMAX && fabs(jet->_etad/10.0)<JETETA)
	{
	  Float_t dRtmp = kinem::delta_R(tau0->_etad,tau0->_phi,jet->_etad/10,jet->_phi);
	  if(dRtmp<dR) {dR = dRtmp; taujet=jet;}
	} 
    } 
    
//     // E N D  R e c a l c u l a t e  M e t  a n d  M e t s i g n i f i c a n c e    
    Int_t Tau_Type = tau0->_tautype;
    Float_t tauEta = tau0->_etad;
    Float_t tauPhi = tau0->_phi;
    Float_t tauPt  = tau0->_Et;
    met_phi = kinem::phi(metx,mety); 
    dphitaumet    = kinem::delta_phi(tau0->_phi, met_phi);
    mttaumet      = TMath::Sqrt(2*tau0->_Et*met*(1-TMath::Cos(dphitaumet)));
    TopTau* tau00 = (TopTau*)taus->At(0);
    if(met>0) ktminp = ktminp*(tau00->_Et+met)/met; else ktminp = -1;

    // F i l l i n g  H i s t o g r a m s F o r QCD/Loose-Tight  s a m p l e
    
    if(sel_taus<=0 && sel_taus_QCD>0 && sel_jetstag>0 && metl>METL){
      weight = 1;
      simu->Fill();
      num_bg++;


    }
  }
  cout<<"num_bg: "<<num_bg<<endl;
  cout << "simu_back->GetEntries() = " << simu->GetEntries() << endl;
  

/////////////////////////////////////BACKGROUND SELECTION ENDS HERE///////////////////////////////////////////////////////////////

//Set1
  TMultiLayerPerceptron *mlp = new TMultiLayerPerceptron("@metl,@ht,@topmassl,@aplan,@sqrts:10:type",
							 "weight",simu,"Entry$%2","Entry$/2");



   mlp->Train(ntrain, "text,graph,update=10");
   // Use TMLPAnalyzer to see what it looks for
   TCanvas* mlpa_canvas = new TCanvas("mlpa_canvas","Network analysis");
   mlpa_canvas->Divide(2,2);
   TMLPAnalyzer* ana = new TMLPAnalyzer(mlp);
   // Initialisation
   ana->GatherInformations();
   // output to the console
   ana->CheckNetwork();
   mlpa_canvas->cd(1);
   // shows how each variable influences the network
   ana->DrawDInputs();
   mlpa_canvas->cd(2);
   // shows the network structure
   mlp->Draw();
   mlpa_canvas->cd(3);
   // draws the resulting network
   ana->DrawNetwork(0,"type==1","type==0");
   mlpa_canvas->cd(4);

   mlpa_canvas->cd(0);
   mlp->Export("Hplus_OUTPUTS/TypeIII/TrainingOutPut/SetI_NNout_SM_type3_tauQCD","C++");
   output->Write();
   output->Close();
   mlpa_canvas->Print("Hplus_OUTPUTS/TypeIII/TrainingOutPut/SetI_NNout_SM_type3_tauQCD.eps");
}