/////////////////
// Markus Klute
// 
// Created: 04/26/2003
// 
// Purpose: Do Mass Analysis
//
#include "top_analyze/D0ChunkAnalyze.hpp"
using namespace TOPexample;
using namespace std;

#include "hitfit/Top_Fit.hpp"
#include "hitfit/Top_Decaykin.hpp"
#include "hitfit/run1_ascii_io.hpp"
#include "hitfit/Lepjets_Event.hpp"
#include "hitfit/Fit_Results.hpp"
#include "hitfit/Fit_Result.hpp"
#include "hitfit/Defaults_Text.hpp"
#include "hitfit/Lepjets_Event_Jet.hpp"
#include "hitfit/Lepjets_Event_Lep.hpp"
#include "hitfit/fourvec.hpp"
#include "CLHEP/Random/JamesRandom.h"

#include <cmath>
#include "kinem_util/AnglesUtil.hpp"
#include "top_tree/TheProperty.hpp"

#include <iostream>

// ============================================================================= 
/*****************************************************************************/
void D0ChunkAnalyze::LJets (std::vector<TheElectronClass>& ElectronVector,
			    std::vector<TheMuonClass>&  MuonVector,
			    std::vector<TheJetClass>& JetVector,
			    float metx, float mety,
			    int run, int evt, float pvz, 
			    bool svxTag, bool sltTag, bool ipTag,
			    int datatype,bool isMC)
{
  // cout<<"use svx was  "<<svxTag<<endl;
  // cout<<"use slt was  "<<sltTag<<endl;
  // cout<<"use ip tag was "<<ipTag<<endl;
 
  std::string df="defaults";
  char* argv[]={"input.dat"}; // does not use it, but has to be declared 
  
  hitfit::Defaults_Text defs (df, 1, argv);
  
  HepJamesRandom engine;
  engine.setSeed (defs.get_int ("seed"));
  
  double lepw_mass  = defs.get_float ("lepw_mass"); 
  double hadw_mass  = defs.get_float ("hadw_mass");
  double top_mass   = defs.get_float ("top_mass");
  int leading_jets  = defs.get_int ("leading_jets");
  double leppt_cut  = defs.get_float ("leppt_cut");
  double jetpt_cut  = defs.get_float ("jetpt_cut");
  double lepeta_cut = defs.get_float ("lepeta_cut");
  double jeteta_cut = defs.get_float ("jeteta_cut");
  double nupt_cut   = defs.get_float ("nupt_cut");
  bool smear        = defs.get_bool ("smear");
  bool print_input_event = defs.get_bool ("print_input_event");
  bool scale_jets   = defs.get_bool ("scale_jets");
  
  // ... instantiate fitter ...	
  hitfit::Top_Fit_Args args (defs);
  hitfit::Top_Fit fitter (args, lepw_mass, hadw_mass, top_mass);
  hitfit::Run1_Ascii_IO_Args io_args (defs);
  
  // ... fill event for fitting ...
  hitfit::Lepjets_Event ev( evt, run);
  ev.setMC(isMC);
  // ...
  // select muon candidate
  std::vector<TheMuonClass>::iterator WMuon;
  int ntmu=0;
  for (std::vector<TheMuonClass>::iterator im = MuonVector.begin(); im != MuonVector.end(); im++)
    {
      if (im->Isolated()|| datatype == 2 )
	{
	  WMuon = im;
	  ntmu++;
	  break;
	}
    }
  
  int nte=0;
  std::vector<TheElectronClass>::iterator WElectron;
  for (std::vector<TheElectronClass>::iterator ie = ElectronVector.begin(); ie != ElectronVector.end(); ie++)
    {
      if (ie->isTight() || datatype == 1 )
	{
	  WElectron = ie;
	  nte++;
	  break;
	}
    }
  
  int lflag = 0;
  if(MuonVector.size()>0      && ElectronVector.size()==0)   lflag=3;
  if(ElectronVector.size()>0  && MuonVector.size() ==0)        lflag=2;
  if( !(ntmu==0 && nte==0) )
    {
      if( ntmu>0 && nte==0 ) lflag=3;
      if( nte>0  && ntmu==0) lflag=2;
      if( (ntmu>0 && nte>0)){
	lflag=2;
	if( WMuon -> pT() > WElectron -> pT()) lflag=3;
      }
    }
  else return;
  
  // ...
  if(lflag==3){ 
    hitfit::Fourvec  lep_fourvec(WMuon -> GetPx(),
				 WMuon -> GetPy(),
				 WMuon -> GetPz(),
				 WMuon -> E());
    ev.add_lep(hitfit::Lepjets_Event_Lep
	       (lep_fourvec,lflag,io_args.muo_res()));
  } 
  
  
  if(lflag==2){
    hitfit::Fourvec  lep_fourvec(WElectron -> GetPx(),
				 WElectron -> GetPy(),
				 WElectron -> GetPz(),
				 WElectron -> E());
    ev.add_lep(hitfit::Lepjets_Event_Lep
	       (lep_fourvec,lflag,io_args.ele_res()));
  } 
  
  float missingET = sqrt ( metx * metx +  mety * mety);
  hitfit::Fourvec met(metx, mety,missingET,missingET);
  ev.met()=met;
  ev.zvertex() = pvz;
  
  std::vector<int> jet_types;
  
  
  bool svx_flag=false;
  int idxjet = -1;
  for ( std::vector<TheJetClass>::iterator ij=JetVector.begin(); ij != JetVector.end() ; ij++)
    {
      idxjet++;

      if(idxjet<4) jet_types.push_back(idxjet+11);
      else jet_types.push_back(20);

      hitfit::Fourvec j_fourvec( ij -> GetPx() , 
				 ij -> GetPy() ,
				 ij -> GetPz() ,
				 ij -> E() );
      bool has_btag = false;
      if ( svxTag && ij->GetSvxTags() == 1 ) has_btag = true;
      if ( sltTag && ij->GetMuoTag() == 1 ) has_btag = true;
      if ( ipTag  && ij->GetIPTag()  == 1 ) has_btag = true;
      ev.add_jet( hitfit::Lepjets_Event_Jet(j_fourvec,jet_types[idxjet],io_args.jet_res(),
					    has_btag,false)  );
    }      
  jet_types . clear();

  ev.kt_res()=io_args.kt_res();
  // ... fill event for fitting done ...
  
  if (smear) ev.smear (engine, true);
  ev.sort ();
  if (print_input_event)  std::cout << ev;
  
  Int_t npermut = 0;
  // ... only events wich pass the cuts ...
  if (ev.cut_leps (leppt_cut, lepeta_cut) < 1 ||
      ev.cut_jets (jetpt_cut, jeteta_cut) < 4 ||
      ev.met().perp() < nupt_cut ||
      ev.nleps() > 1)
    {
      std::cout << "failed kin cuts\n";
    }
  else {
    
    if (leading_jets > 0) {
      assert (leading_jets >= 4);
      ev.trimjets (leading_jets);
    }
    
    hitfit::Fit_Results fitres = fitter.fit (ev);
    if (fitres[hitfit::all_list].size() > 0) {
      
      //std::cout << fitres[hitfit::all_list][0].chisq()
      //   << " " << fitres[hitfit::all_list][0].mt()
      //   << " " << fitres[hitfit::all_list][0].sigmt()
      //   << "\n";
      
      npermut =13;
      Float_t  chisq[13];  // Fit chisq.< 0 means fit didn't converge.
      Float_t  mt[13];     // Top mass from fit.
      Float_t  sigmt[13];  // Top mass uncertainty from fit
      Float_t  umwhad[13]; // Mass of hadronic W in this permutation, before fitting.
      Float_t  utmass[13]; // Mass of hadronic top in this permutation, before fitting.
      Int_t    perm[13];   // Jet types, packed 4 bits per jets.
      Float_t  m_tt[13];   // Invariant mass of ttbar pair.
      Float_t  pt_t[13];   // Top pt, averaged over both.
      Float_t  kt[13];     // Event kT.
      Float_t  eta_t1[13]; // Eta of leptonic top.
      Float_t  eta_t2[13]; // Eta of hadronic top.
      Float_t  phi_t1[13]; // Phi of leptonic top.
      Float_t  phi_t2[13]; // Phi of hadronic top.
      Float_t  pt_t1[13];  // Pt of leptonic top.
      Float_t  pt_t2[13];  // Pt of hadronic top.

      Float_t t1[13][3];   // eta,phi,pt of leptonic top
      Float_t t2[13][3];   // eta,phi,pt of hadronic top

      Float_t w1[13][3];   // eta,phi,pt of leptonic W
      Float_t w2[13][3];   // eta,phi,pt of hadronic W

      Float_t jet_lepb[13][4];  //px,py,pz,E of b-jet from leptonic Top after fitting
      Float_t jet_hadb[13][4];  //px,py,pz,E of b-jet from hadronic Top after fitting
      Float_t jet_hadw1[13][4]; //px,py,pz,E of first jet from hadronic W after fitting
      Float_t jet_hadw2[13][4]; //px,py,pz,E of second jet from hadronic W after fitting
      Float_t lepton[13][4];    //px,py,pz,E of Lepton after fitting
      Float_t neutrino[13][4];  //px,py,pz,E of Neutrino after fitting

      //the same but with btagging
      Float_t  bchisq[13];  // Fit chisq.< 0 means fit didn't converge.
      Float_t  bmt[13];     // Top mass from fit.
      Float_t  bsigmt[13];  // Top mass uncertainty from fit
      Float_t  bumwhad[13]; // Mass of hadronic W in this permutation, before fitting.
      Float_t  butmass[13]; // Mass of hadronic top in this permutation, before fitting.
      Int_t    bperm[13];   // Jet types, packed 4 bits per jets.
      Float_t  bm_tt[13];   // Invariant mass of ttbar pair.
      Float_t  bpt_t[13];   // Top pt, averaged over both.
      Float_t  bkt[13];     // Event kT.
      Float_t  beta_t1[13]; // Eta of leptonic top.
      Float_t  beta_t2[13]; // Eta of hadronic top.
      Float_t  bphi_t1[13]; // Phi of leptonic top.
      Float_t  bphi_t2[13]; // Phi of hadronic top.
      Float_t  bpt_t1[13];  // Pt of leptonic top.
      Float_t  bpt_t2[13];  // Pt of hadronic top.

      Float_t bt1[13][3];   // eta,phi,pt of leptonic top
      Float_t bt2[13][3];   // eta,phi,pt of hadronic top

      Float_t bw1[13][3];   // eta,phi,pt of leptonic W
      Float_t bw2[13][3];   // eta,phi,pt of hadronic W

      Float_t bjet_lepb[13][4];  //px,py,pz,E of b-jet from leptonic Top after fitting
      Float_t bjet_hadb[13][4];  //px,py,pz,E of b-jet from hadronic Top after fitting
      Float_t bjet_hadw1[13][4]; //px,py,pz,E of first jet from hadronic W after fitting
      Float_t bjet_hadw2[13][4]; //px,py,pz,E of second jet from hadronic W after fitting
      Float_t blepton[13][4];    //px,py,pz,E of Lepton after fitting
      Float_t bneutrino[13][4];  //px,py,pz,E of Neutrino after fitting
      
      for(int i=0; i<npermut; i++){
	chisq[i] = -10.;
	mt[i] = sigmt[i] = umwhad[i] = utmass[i] = perm[i] = 0.; 
	m_tt[i] = pt_t[i] = kt[i] = 0.;
	for(int j=0;j<3;j++){
	  t1[i][j]  = t2[i][j]  = w1[i][j]  = w2[i][j]  = 0.;
	  bt1[i][j] = bt2[i][j] = bw1[i][j] = bw2[i][j] = 0.;
	}
	for(int j=0;j<4;j++){
	  jet_lepb[i][j]  = jet_hadb[i][j]  = jet_hadw1[i][j]  = jet_hadw2[i][j]=0.;
	  lepton[i][j]=neutrino[i][j]=0.;
	  bjet_lepb[i][j] = bjet_hadb[i][j] = bjet_hadw1[i][j] = bjet_hadw2[i][j]=0.;
	  blepton[i][j]=bneutrino[i][j]=0.;
	}
	bchisq[i] = -10.;
	bmt[i] = bsigmt[i] = bumwhad[i] = butmass[i] = bperm[i] = 0.; 
	bm_tt[i] = bpt_t[i] = bkt[i]  = 0.;
      }
      
      // ... get the fit variables ...
      int packed = 0;
      int bpacked = 0;
      if (fitres[hitfit::noperm_list].size() > 0) 
	{
	  chisq[0] =fitres[hitfit::noperm_list][0].chisq();
	  mt[0]    =fitres[hitfit::noperm_list][0].mt();
	  sigmt[0] =fitres[hitfit::noperm_list][0].sigmt();
	  umwhad[0]=fitres[hitfit::noperm_list][0].umwhad();
	  utmass[0]=fitres[hitfit::noperm_list][0].utmass();
	  const hitfit::Lepjets_Event& ev_fitted = fitres[hitfit::noperm_list][0].ev();
	  packed = 0;
	  for (int ii=0; ii < ev.njets(); ii++)
	    packed = ((packed << 4) | (ev_fitted.jet(ii).type() & 0xf));
	  
	  perm[0]=packed;
	  hitfit::Fourvec hadt = hitfit::Top_Decaykin::hadt(ev_fitted);
	  hitfit::Fourvec lept = hitfit::Top_Decaykin::lept(ev_fitted);	      
	  hitfit::Fourvec hadw = hitfit::Top_Decaykin::hadw(ev_fitted);
	  hitfit::Fourvec lepw = hitfit::Top_Decaykin::lepw(ev_fitted);		  
	  
	  m_tt[0]   = (hadt + lept).m ();
	  pt_t[0]   = (hadt.perp() + lept.perp()) / 2; 
	  kt[0]     = ev_fitted.kt().perp();   
	  t1[0][0]  = lept.pseudoRapidity();
	  t2[0][0]  = hadt.pseudoRapidity();
	  t1[0][1]  = lept.phi();
	  t2[0][1]  = hadt.phi();
	  t1[0][2]  = lept.perp(); 
	  t2[0][2]  = hadt.perp(); 
	  w1[0][0]  = lepw.pseudoRapidity();
	  w2[0][0]  = hadw.pseudoRapidity();
	  w1[0][1]  = lepw.phi();
	  w2[0][1]  = hadw.phi();
	  w1[0][2]  = lepw.perp(); 
	  w2[0][2]  = hadw.perp(); 



	  for (int ii=0; ii < ev.njets(); ii++){
	    hitfit::Fourvec rjet=ev_fitted.jet(ii).p();
	    switch(ev_fitted.jet(ii).type()){
	      
	    case hitfit::lepb_label : 
	      jet_lepb[0][0]=rjet.px();
	      jet_lepb[0][1]=rjet.py();
	      jet_lepb[0][2]=rjet.pz();
	      jet_lepb[0][3]=rjet.e();
	      break;
	      
	    case hitfit::hadb_label : 
	      jet_hadb[0][0]=rjet.px();
	      jet_hadb[0][1]=rjet.py();
	      jet_hadb[0][2]=rjet.pz();
	      jet_hadb[0][3]=rjet.e();
	      break;
	      
	    case hitfit::hadw1_label : 
	      jet_hadw1[0][0]=rjet.px();
	      jet_hadw1[0][1]=rjet.py();
	      jet_hadw1[0][2]=rjet.pz();
	      jet_hadw1[0][3]=rjet.e();
	      break;
	      
	    case hitfit::hadw2_label : 
	      jet_hadw2[0][0]=rjet.px();
	      jet_hadw2[0][1]=rjet.py();
	      jet_hadw2[0][2]=rjet.pz();
	      jet_hadw2[0][3]=rjet.e();
	      break;
	      
	    default:
	      //  cerr<<"Unknown Jet-Label in HitFit: "<<ev_fitted_i.jet(ii).type()<<endl;
	      break;
	    }
	  }

	  hitfit::Fourvec rlep=ev_fitted.lep(0).p();
	  hitfit::Fourvec rnu =ev_fitted.met();

	  lepton[0][0]=rlep.px();
	  lepton[0][1]=rlep.py();
	  lepton[0][2]=rlep.pz();
	  lepton[0][3]=rlep.e();

	  neutrino[0][0]=rnu.px();
	  neutrino[0][1]=rnu.py();
	  neutrino[0][2]=rnu.pz();
	  neutrino[0][3]=rnu.e();
	  



	  // cout<<"FINISHED GETTING 'correct' FIT INFO WITHOUT BTAGGING"<<endl;
	}
            
      int solutions = fitres[hitfit::all_list].size();
      
      for(int i=1;(i < npermut && i <solutions + 1); i++)
	{
	  chisq[i] =fitres[hitfit::all_list][i-1].chisq();
	  mt[i]    =fitres[hitfit::all_list][i-1].mt();
	  sigmt[i] =fitres[hitfit::all_list][i-1].sigmt();
	  umwhad[i]=fitres[hitfit::all_list][i-1].umwhad();
	  utmass[i]=fitres[hitfit::all_list][i-1].utmass();
	  const hitfit::Lepjets_Event& ev_fitted_i = fitres[hitfit::all_list][i-1].ev();
	  packed = 0;
	  for (int ii=0; ii < ev.njets(); ii++)
	    packed = ((packed << 4) | (ev_fitted_i.jet(ii).type() & 0xf));
	  perm[i]=packed;
	  
	  hitfit::Fourvec hadt = hitfit::Top_Decaykin::hadt(ev_fitted_i);
	  hitfit::Fourvec lept = hitfit::Top_Decaykin::lept(ev_fitted_i);
	      
	  hitfit::Fourvec hadw = hitfit::Top_Decaykin::hadw(ev_fitted_i);
	  hitfit::Fourvec lepw = hitfit::Top_Decaykin::lepw(ev_fitted_i);
	      
	  m_tt[i]   = (hadt + lept).m ();
	  pt_t[i]   = (hadt.perp() + lept.perp()) / 2; 
	  kt[i]     = ev_fitted_i.kt().perp();   

	  t1[i][0]  = lept.pseudoRapidity();
	  t2[i][0]  = hadt.pseudoRapidity();
	  t1[i][1]  = lept.phi();
	  t2[i][1]  = hadt.phi();
	  t1[i][2]  = lept.perp(); 
	  t2[i][2]  = hadt.perp(); 

	  w1[i][0]  = lepw.pseudoRapidity();
	  w2[i][0]  = hadw.pseudoRapidity();
	  w1[i][1]  = lepw.phi();
	  w2[i][1]  = hadw.phi();
	  w1[i][2]  = lepw.perp(); 
	  w2[i][2]  = hadw.perp(); 

	  for (int ii=0; ii < ev.njets(); ii++){
	    hitfit::Fourvec rjet=ev_fitted_i.jet(ii).p();
	    switch(ev_fitted_i.jet(ii).type()){
	      
	    case hitfit::lepb_label : 
	      jet_lepb[i][0]=rjet.px();
	      jet_lepb[i][1]=rjet.py();
	      jet_lepb[i][2]=rjet.pz();
	      jet_lepb[i][3]=rjet.e();
	      break;
	      
	    case hitfit::hadb_label : 
	      jet_hadb[i][0]=rjet.px();
	      jet_hadb[i][1]=rjet.py();
	      jet_hadb[i][2]=rjet.pz();
	      jet_hadb[i][3]=rjet.e();
	      break;
	      
	    case hitfit::hadw1_label : 
	      jet_hadw1[i][0]=rjet.px();
	      jet_hadw1[i][1]=rjet.py();
	      jet_hadw1[i][2]=rjet.pz();
	      jet_hadw1[i][3]=rjet.e();
	      break;
	      
	    case hitfit::hadw2_label : 
	      jet_hadw2[i][0]=rjet.px();
	      jet_hadw2[i][1]=rjet.py();
	      jet_hadw2[i][2]=rjet.pz();
	      jet_hadw2[i][3]=rjet.e();
	      break;
	      
	    default:
	      //  cerr<<"Unknown Jet-Label in HitFit: "<<ev_fitted_i.jet(ii).type()<<endl;
	      break;
	    }
	  }

	  hitfit::Fourvec rlep=ev_fitted_i.lep(0).p();
	  hitfit::Fourvec rnu =ev_fitted_i.met();

	  lepton[i][0]=rlep.px();
	  lepton[i][1]=rlep.py();
	  lepton[i][2]=rlep.pz();
	  lepton[i][3]=rlep.e();

	  neutrino[i][0]=rnu.px();
	  neutrino[i][1]=rnu.py();
	  neutrino[i][2]=rnu.pz();
	  neutrino[i][3]=rnu.e();
	  
	  
	}
      // cout<<"FINISHED GETTING ALL FITS INFO WITHOUT BTAGGING"<<endl;
      int bsolutions = fitres[hitfit::btag_list].size();
      // cout<<"For this event there were " <<bsolutions<<"  solutions consistent with b tags"<<endl;
      for(int i = 1; (i < npermut && i < bsolutions +1); i++) {
	bchisq[i] =fitres[hitfit::btag_list][i-1].chisq();
	bmt[i]    =fitres[hitfit::btag_list][i-1].mt();
	bsigmt[i] =fitres[hitfit::btag_list][i-1].sigmt();
	bumwhad[i]=fitres[hitfit::btag_list][i-1].umwhad();
	butmass[i]=fitres[hitfit::btag_list][i-1].utmass();
	const hitfit::Lepjets_Event& ev_fitted_ib = fitres[hitfit::btag_list][i-1].ev();
	bpacked = 0;
	for (int ii=0; ii < ev.njets(); ii++)
	  bpacked = ((bpacked << 4) | (ev_fitted_ib.jet(ii).type() & 0xf));
	bperm[i]=bpacked;
	
	hitfit::Fourvec hadtb = hitfit::Top_Decaykin::hadt(ev_fitted_ib);
	hitfit::Fourvec leptb = hitfit::Top_Decaykin::lept(ev_fitted_ib);	      

	hitfit::Fourvec hadwb = hitfit::Top_Decaykin::hadw(ev_fitted_ib);
	hitfit::Fourvec lepwb = hitfit::Top_Decaykin::lepw(ev_fitted_ib);

	bm_tt[i]   = (hadtb + leptb).m ();
	bpt_t[i]   = (hadtb.perp() + leptb.perp()) / 2; 
	bkt[i]     = ev_fitted_ib.kt().perp();   

	bt1[i][0] = leptb.pseudoRapidity();
	bt2[i][0] = hadtb.pseudoRapidity();
	bt1[i][1] = leptb.phi();
	bt2[i][1] = hadtb.phi();
	bt1[i][2]  = leptb.perp(); 
	bt2[i][2]  = hadtb.perp(); 

	bw1[i][0] = lepwb.pseudoRapidity();
	bw2[i][0] = hadwb.pseudoRapidity();
	bw1[i][1] = lepwb.phi();
	bw2[i][1] = hadwb.phi();
	bw1[i][2]  = lepwb.perp(); 
	bw2[i][2]  = hadwb.perp(); 

	for (int ii=0; ii < ev.njets(); ii++){
	  hitfit::Fourvec rjet=ev_fitted_ib.jet(ii).p();
	  switch(ev_fitted_ib.jet(ii).type()){
	      
	    case hitfit::lepb_label : 
	      bjet_lepb[i][0]=rjet.px();
	      bjet_lepb[i][1]=rjet.py();
	      bjet_lepb[i][2]=rjet.pz();
	      bjet_lepb[i][3]=rjet.e();
	      break;
	      
	    case hitfit::hadb_label : 
	      bjet_hadb[i][0]=rjet.px();
	      bjet_hadb[i][1]=rjet.py();
	      bjet_hadb[i][2]=rjet.pz();
	      bjet_hadb[i][3]=rjet.e();
	      break;
	      
	    case hitfit::hadw1_label : 
	      bjet_hadw1[i][0]=rjet.px();
	      bjet_hadw1[i][1]=rjet.py();
	      bjet_hadw1[i][2]=rjet.pz();
	      bjet_hadw1[i][3]=rjet.e();
	      break;
	      
	    case hitfit::hadw2_label : 
	      bjet_hadw2[i][0]=rjet.px();
	      bjet_hadw2[i][1]=rjet.py();
	      bjet_hadw2[i][2]=rjet.pz();
	      bjet_hadw2[i][3]=rjet.e();
	      break;
	      
	    default:
	      //  cerr<<"Unknown Jet-Label in HitFit: "<<ev_fitted_i.jet(ii).type()<<endl;
	      break;
	    }
	  }

	  hitfit::Fourvec rlep=ev_fitted_ib.lep(0).p();
	  hitfit::Fourvec rnu =ev_fitted_ib.met();

	  blepton[i][0]=rlep.px();
	  blepton[i][1]=rlep.py();
	  blepton[i][2]=rlep.pz();
	  blepton[i][3]=rlep.e();

	  bneutrino[i][0]=rnu.px();
	  bneutrino[i][1]=rnu.py();
	  bneutrino[i][2]=rnu.pz();
	  bneutrino[i][3]=rnu.e();

	
      }
      // cout<<"FINISHED GETTING ALL FITS INFO WITH BTAGGING"<<endl;
      for(int i=0;(i < npermut && i <solutions + 1); i++){
	_theProperty -> AddTopFit
	  (chisq[i], mt[i], sigmt[i], 
	   umwhad[i], utmass[i], m_tt[i], pt_t[i], kt[i], 
	   t1[i], t2[i], perm[i], w1[i], w2[i],
	   jet_lepb[i],jet_hadb[i],jet_hadw1[i], jet_hadw2[i],
	   lepton[i], neutrino[i],
	   bchisq[i], bmt[i], bsigmt[i], 
	   bumwhad[i], butmass[i], bm_tt[i], bpt_t[i], bkt[i], 
	   bt1[i], bt2[i], bperm[i], bw1[i], bw2[i],
	   bjet_lepb[i], bjet_hadb[i], bjet_hadw1[i], bjet_hadw2[i],
	   blepton[i], bneutrino[i]
	   );
      }
    } 
  }
 
}
/*****************************************************************************/