#include "makeTree.h"

makeTree::makeTree(TString sample)
{

  if (sample.Contains("data"))
    isRealData=true;
  else
    isRealData=false;

  sample_ = sample;

  cout<<"constructing tree for sample: "<<sample<<endl; 

  //input to copy/create branches from
  TString filename = "/cu3/galank/DileptonTopReconstruction/" + sample + "/edmntuple.root";
  fin = TFile::Open(filename);

  //output
  filename = "/cu3/galank/DileptonTopReconstruction/reducedTrees/reducedTree_" + sample + ".root";
  fout = new TFile(filename,"RECREATE","reducedTree");
  fout->cd();

  tree = new TTree("tree","tree");

  tree->Branch( "lightJet1Pt"           ,  &lightJet1Pt           ); 
  tree->Branch( "lightJet1Eta"          ,  &lightJet1Eta          );
  tree->Branch( "lightJet1Phi"          ,  &lightJet1Phi          );
  tree->Branch( "lightJet2Pt"           ,  &lightJet2Pt           ); 
  tree->Branch( "lightJet2Eta"          ,  &lightJet2Eta          );
  tree->Branch( "lightJet2Phi"          ,  &lightJet2Phi          );
  tree->Branch( "bJet1Pt"               ,  &bJet1Pt               );
  tree->Branch( "bJet1Eta"              ,  &bJet1Eta              );
  tree->Branch( "bJet1Phi"              ,  &bJet1Phi              );
  tree->Branch( "bJet2Pt"               ,  &bJet2Pt               ); 
  tree->Branch( "bJet2Eta"              ,  &bJet2Eta              );
  tree->Branch( "bJet2Phi"              ,  &bJet2Phi              );
  tree->Branch( "nLoosebJets"           ,  &nLoosebJets           );
  tree->Branch( "nTightbJets"           ,  &nTightbJets           );
  tree->Branch( "lepton1Pt"             ,  &lepton1Pt             );
  tree->Branch( "lepton1Eta"            ,  &lepton1Eta            );
  tree->Branch( "lepton1Phi"            ,  &lepton1Phi            );
  tree->Branch( "lepton1Flavor"         ,  &lepton1Flavor         );
  tree->Branch( "lepton1Charge"         ,  &lepton1Charge         );
  tree->Branch( "lepton2Pt"             ,  &lepton2Pt             ); 
  tree->Branch( "lepton2Eta"            ,  &lepton2Eta            );
  tree->Branch( "lepton2Phi"            ,  &lepton2Phi            );
  tree->Branch( "lepton2Flavor"         ,  &lepton2Flavor         );
  tree->Branch( "lepton2Charge"         ,  &lepton2Charge         );
  tree->Branch(  "doubleEle"            ,  &doubleEle             );
  tree->Branch(  "doubleMu"             ,  &doubleMu              );
  tree->Branch( "sameFlavorLeptons"     ,  &sameFlavorLeptons     );
  tree->Branch( "inZwindow"             ,  &inZwindow             );
  tree->Branch( "dileptonInvariantMass" ,  &dileptonInvariantMass );
  tree->Branch( "nElectrons"            ,  &nElectrons            );
  tree->Branch( "nMuons"                ,  &nMuons                );
  tree->Branch( "muon1Pt"               ,  &muon1Pt               );
  tree->Branch( "muon1Eta"              ,  &muon1Eta              );
  tree->Branch( "ele1Pt"                ,  &ele1Pt                );
  tree->Branch( "ele1Eta"               ,  &ele1Eta               );
  tree->Branch( "MET"                   ,  &MET                   );
  tree->Branch( "kinCut"                ,  &kinCut                );
  tree->Branch( "Zveto"                 ,  &Zveto                 );
  tree->Branch( "wideInvertZveto"       ,  &wideInvertZveto       );
  tree->Branch( "invertZveto"           ,  &invertZveto           );

  //MuHad triggers
  tree->Branch( "HLT_IsoMu17_eta2p1_TriCentralPFJet30"                   ,  &HLT_IsoMu17_eta2p1_TriCentralPFJet30                   );	//		 
  tree->Branch( "HLT_IsoMu20_eta2p1_CentralPFJet30_BTagIPIter"           ,  &HLT_IsoMu20_eta2p1_CentralPFJet30_BTagIPIter           );		 
  tree->Branch( "HLT_IsoMu20_eta2p1_CentralPFJet80"                      ,  &HLT_IsoMu20_eta2p1_CentralPFJet80                      );	//		 
  tree->Branch( "HLT_IsoMu20_eta2p1_CentralPFNoPUJet30_BTagIPIter"       ,  &HLT_IsoMu20_eta2p1_CentralPFNoPUJet30_BTagIPIter       );	 
  tree->Branch( "HLT_IsoMu20_eta2p1_TriCentralPFJet30"                   ,  &HLT_IsoMu20_eta2p1_TriCentralPFJet30                   );	//		 
  tree->Branch( "HLT_IsoMu20_eta2p1_TriCentralPFNoPUJet30"               ,  &HLT_IsoMu20_eta2p1_TriCentralPFNoPUJet30               );		 
  tree->Branch( "HLT_IsoMu24_eta2p1_CentralPFJet30_CentralPFJet25"       ,  &HLT_IsoMu24_eta2p1_CentralPFJet30_CentralPFJet25       );	 //
  tree->Branch( "HLT_Mu15_eta2p1_TriCentral_40_20_20_DiBTagIP3D1stTrack" ,  &HLT_Mu15_eta2p1_TriCentral_40_20_20_DiBTagIP3D1stTrack );
  tree->Branch( "HLT_Mu8_Ele8_CaloIdT_TrkIdVL_Mass8_PFHT175"             ,  &HLT_Mu8_Ele8_CaloIdT_TrkIdVL_Mass8_PFHT175             );		 
  tree->Branch( "HLT_Mu8_TriJet30"                                       ,  &HLT_Mu8_TriJet30                                       );
  //DoubleMu triggers
  tree->Branch( "HLT_DoubleMu5_IsoMu5",  &HLT_DoubleMu5_IsoMu5 );
  tree->Branch( "HLT_Mu17_Mu8"        ,  &HLT_Mu17_Mu8         );
  tree->Branch( "HLT_Mu17_TkMu8"      ,  &HLT_Mu17_TkMu8       );
  tree->Branch( "HLT_Mu22_TkMu8"      ,  &HLT_Mu22_TkMu8       );
  tree->Branch( "HLT_Mu22_TkMu22"     ,  &HLT_Mu22_TkMu22      );
  //DoubleElectron triggers
  tree->Branch( "HLT_DoubleEle10_CaloIdL_TrkIdVL_Ele10_CaloIdT_TrkIdVL" ,  &HLT_DoubleEle10_CaloIdL_TrkIdVL_Ele10_CaloIdT_TrkIdVL );
  tree->Branch( "HLT_Ele15_Ele8_Ele5_CaloIdL_TrkIdVL"                   ,  &HLT_Ele15_Ele8_Ele5_CaloIdL_TrkIdVL                   );
  tree->Branch( "HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL", &HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL );
  //MuEG triggers
  tree->Branch( "HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL",  &HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL );
  tree->Branch( "HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL",  &HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL );


  sbottomTree = new TTree("sbottomTree","sbottomTree");
  sbottomTree->Branch( "foundPair"             ,  &foundPair       );
  sbottomTree->Branch( "sbottomCandMass"       ,  &sbottomCandMass );
  sbottomTree->Branch( "sbottomTruth"          ,  &sbottomTruth    );
  sbottomTree->Branch( "deltaM"                ,  &deltaM          );
  sbottomTree->Branch( "lightJet1Pt"           ,  &lightJet1Pt           );
  sbottomTree->Branch( "lightJet1Eta"          ,  &lightJet1Eta          );
  sbottomTree->Branch( "lightJet2Pt"           ,  &lightJet2Pt           );
  sbottomTree->Branch( "lightJet2Eta"          ,  &lightJet2Eta          );
  sbottomTree->Branch( "bJet1Pt"               ,  &bJet1Pt               );
  sbottomTree->Branch( "bJet1Eta"              ,  &bJet1Eta              );
  sbottomTree->Branch( "bJet2Pt"               ,  &bJet2Pt               );
  sbottomTree->Branch( "bJet2Eta"              ,  &bJet2Eta              );
  sbottomTree->Branch( "lepton1Pt"             ,  &lepton1Pt             );
  sbottomTree->Branch( "lepton1Eta"            ,  &lepton1Eta            );
  sbottomTree->Branch( "lepton2Pt"             ,  &lepton2Pt             );
  sbottomTree->Branch( "lepton2Eta"            ,  &lepton2Eta            );
  sbottomTree->Branch( "sameFlavorLeptons"     ,  &sameFlavorLeptons     );
  sbottomTree->Branch( "inZwindow"             ,  &inZwindow             );
  sbottomTree->Branch( "dileptonInvariantMass" ,  &dileptonInvariantMass );
}

makeTree::~makeTree(){

}

void makeTree::fillTree()
{

  if (fin) {
    fwlite::Event ev(fin);
    cout<<"Beginning loop over events..."<<endl;

    int nEv=0;

    for (ev.toBegin(); !ev.atEnd(); ++ev) {

      //cout<<ev.eventAuxiliary().run()<<" "<<ev.eventAuxiliary().luminosityBlock()<<" "<<ev.eventAuxiliary().event()<<endl;

      nEv++;
      //if (nEv>100) break;
      if (nEv % 10000 == 0) cout<<"filling event: "<<nEv<<endl;

      //for data
      fwlite::Handle< int > nlightJets;
      nlightJets.getByLabel(ev,"edmNtuple","nLightJets");
      if (!nlightJets.isValid()) {
        cout<<"invalid handle nLightJets in event: "<<ev.eventAuxiliary().run()<<" "<<ev.eventAuxiliary().luminosityBlock()<<" "<<ev.eventAuxiliary().event()<<endl;
        continue;
      }


      //trigger
      if (isRealData) {
	HLT_IsoMu17_eta2p1_TriCentralPFJet30=false;
	HLT_IsoMu20_eta2p1_CentralPFJet30_BTagIPIter=false;
	HLT_IsoMu20_eta2p1_CentralPFJet80=false;
        HLT_IsoMu20_eta2p1_CentralPFNoPUJet30_BTagIPIter=false;
        HLT_IsoMu20_eta2p1_TriCentralPFJet30=false;
        HLT_IsoMu20_eta2p1_TriCentralPFNoPUJet30=false;
        HLT_IsoMu24_eta2p1_CentralPFJet30_CentralPFJet25=false;
        HLT_Mu15_eta2p1_TriCentral_40_20_20_DiBTagIP3D1stTrack=false;
        HLT_Mu8_Ele8_CaloIdT_TrkIdVL_Mass8_PFHT175=false;
        HLT_Mu8_TriJet30=false;

	HLT_DoubleMu5_IsoMu5=false;
	HLT_Mu17_Mu8=false;
	HLT_Mu17_TkMu8=false;
	HLT_Mu22_TkMu8=false;
	HLT_Mu22_TkMu22=false;

	HLT_DoubleEle10_CaloIdL_TrkIdVL_Ele10_CaloIdT_TrkIdVL=false;
	HLT_Ele15_Ele8_Ele5_CaloIdL_TrkIdVL=false;
	HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL=false;

	HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL=false;
	HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL=false;

	//cout<<"adding trigger information"<<endl;

	fwlite::Handle< edm::TriggerResults > trigRes;
        trigRes.getByLabel(ev,"TriggerResults","","HLT");
	edm::TriggerResults triggerResults = trigRes.ref();
	edm::TriggerNames triggerNames = ev.triggerNames(triggerResults);
        for (unsigned int i=0; i<triggerNames.size(); i++) {
          string name = triggerNames.triggerName(i);
	  if ((TString(name)).Contains( "HLT_IsoMu17_eta2p1_TriCentralPFJet30_v" ))
	    HLT_IsoMu17_eta2p1_TriCentralPFJet30 = triggerResults.accept(triggerNames.triggerIndex(name));			
	  if ((TString(name)).Contains( "HLT_IsoMu20_eta2p1_CentralPFJet30_BTagIPIter_v" ))
	    HLT_IsoMu20_eta2p1_CentralPFJet30_BTagIPIter = triggerResults.accept(triggerNames.triggerIndex(name));		
	  if ((TString(name)).Contains( "HLT_IsoMu20_eta2p1_CentralPFJet80_v" ))
	    HLT_IsoMu20_eta2p1_CentralPFJet80 = triggerResults.accept(triggerNames.triggerIndex(name));			
	  if ((TString(name)).Contains( "HLT_IsoMu20_eta2p1_CentralPFNoPUJet30_BTagIPIter_v" ))
	    HLT_IsoMu20_eta2p1_CentralPFNoPUJet30_BTagIPIter = triggerResults.accept(triggerNames.triggerIndex(name));	
	  if ((TString(name)).Contains( "HLT_IsoMu20_eta2p1_TriCentralPFJet30_v" ))
	    HLT_IsoMu20_eta2p1_TriCentralPFJet30 = triggerResults.accept(triggerNames.triggerIndex(name));			
	  if ((TString(name)).Contains( "HLT_IsoMu20_eta2p1_TriCentralPFNoPUJet30_v" ))
	    HLT_IsoMu20_eta2p1_TriCentralPFNoPUJet30 = triggerResults.accept(triggerNames.triggerIndex(name));		
	  if ((TString(name)).Contains( "HLT_IsoMu24_eta2p1_CentralPFJet30_CentralPFJet25_v" ))
	    HLT_IsoMu24_eta2p1_CentralPFJet30_CentralPFJet25 = triggerResults.accept(triggerNames.triggerIndex(name));	
	  if ((TString(name)).Contains( "HLT_Mu15_eta2p1_TriCentral_40_20_20_DiBTagIP3D1stTrack_v" ))
	    HLT_Mu15_eta2p1_TriCentral_40_20_20_DiBTagIP3D1stTrack = triggerResults.accept(triggerNames.triggerIndex(name));
	  if ((TString(name)).Contains( "HLT_Mu8_Ele8_CaloIdT_TrkIdVL_Mass8_PFHT175_v" ))
	    HLT_Mu8_Ele8_CaloIdT_TrkIdVL_Mass8_PFHT175 = triggerResults.accept(triggerNames.triggerIndex(name));		
	  if ((TString(name)).Contains( "HLT_Mu8_TriJet30_v" ))
	    HLT_Mu8_TriJet30 = triggerResults.accept(triggerNames.triggerIndex(name));                                      

	  if ((TString(name)).Contains( "HLT_DoubleMu5_IsoMu5" ))
	    HLT_DoubleMu5_IsoMu5 = triggerResults.accept(triggerNames.triggerIndex(name));
	  if ((TString(name)).Contains( "HLT_Mu17_Mu8" ))
	    HLT_Mu17_Mu8 = triggerResults.accept(triggerNames.triggerIndex(name));
	  if ((TString(name)).Contains( "HLT_Mu17_TkMu8" ))
            HLT_Mu17_TkMu8 = triggerResults.accept(triggerNames.triggerIndex(name));
	  if ((TString(name)).Contains( "HLT_Mu22_TkMu8" ))
	    HLT_Mu22_TkMu8 = triggerResults.accept(triggerNames.triggerIndex(name));
	  if ((TString(name)).Contains( "HLT_Mu22_TkMu22" ))
	    HLT_Mu22_TkMu22 = triggerResults.accept(triggerNames.triggerIndex(name));
	  
	  if ((TString(name)).Contains( "HLT_DoubleEle10_CaloIdL_TrkIdVL_Ele10_CaloIdT_TrkIdVL" ))
	    HLT_DoubleEle10_CaloIdL_TrkIdVL_Ele10_CaloIdT_TrkIdVL = triggerResults.accept(triggerNames.triggerIndex(name));
	  if ((TString(name)).Contains( "HLT_Ele15_Ele8_Ele5_CaloIdL_TrkIdVL" ))
	    HLT_Ele15_Ele8_Ele5_CaloIdL_TrkIdVL = triggerResults.accept(triggerNames.triggerIndex(name));
	  if ((TString(name)).Contains( "HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL" ))
	    HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL = triggerResults.accept(triggerNames.triggerIndex(name));
	  
	  if ((TString(name)).Contains( "HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL" ))
	    HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL = triggerResults.accept(triggerNames.triggerIndex(name));
	  if ((TString(name)).Contains( "HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL" ))
	    HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL = triggerResults.accept(triggerNames.triggerIndex(name));

	} //loop over triggers
      } //data
      else {
	HLT_IsoMu17_eta2p1_TriCentralPFJet30=true;
	HLT_IsoMu20_eta2p1_CentralPFJet30_BTagIPIter=true;
	HLT_IsoMu20_eta2p1_CentralPFJet80=true;
	HLT_IsoMu20_eta2p1_CentralPFNoPUJet30_BTagIPIter=true;
	HLT_IsoMu20_eta2p1_TriCentralPFJet30=true;
	HLT_IsoMu20_eta2p1_TriCentralPFNoPUJet30=true;
	HLT_IsoMu24_eta2p1_CentralPFJet30_CentralPFJet25=true;
	HLT_Mu15_eta2p1_TriCentral_40_20_20_DiBTagIP3D1stTrack=true;
	HLT_Mu8_Ele8_CaloIdT_TrkIdVL_Mass8_PFHT175=true;
	HLT_Mu8_TriJet30=true;

	HLT_DoubleMu5_IsoMu5=true;
        HLT_Mu17_Mu8=true;
	HLT_Mu17_TkMu8=true;
        HLT_Mu22_TkMu8=true;
	HLT_Mu22_TkMu22=true;

        HLT_DoubleEle10_CaloIdL_TrkIdVL_Ele10_CaloIdT_TrkIdVL=true;
        HLT_Ele15_Ele8_Ele5_CaloIdL_TrkIdVL=true;
	HLT_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL=true;

        HLT_Mu17_Ele8_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL=true;
        HLT_Mu8_Ele17_CaloIdT_CaloIsoVL_TrkIdVL_TrkIsoVL=true;
      } //MC


      //cout<<"light jets"<<endl;


      //light jets
      fwlite::Handle< ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > lightJet1Vector;
      lightJet1Vector.getByLabel(ev,"edmNtuple","lightJet1LorentzVector");
      lightJet1.SetPxPyPzE(lightJet1Vector->Px(),lightJet1Vector->Py(),lightJet1Vector->Pz(),lightJet1Vector->E());
      fwlite::Handle< ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > lightJet2Vector;
      lightJet2Vector.getByLabel(ev,"edmNtuple","lightJet2LorentzVector");
      lightJet2.SetPxPyPzE(lightJet2Vector->Px(),lightJet2Vector->Py(),lightJet2Vector->Pz(),lightJet2Vector->E());


      //cout<<"leptons"<<endl;


      //leptons
      fwlite::Handle< ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > lepton1Vector;
      lepton1Vector.getByLabel(ev,"edmNtuple","lepton1LorentzVector");
      TLorentzVector lepton1; lepton1.SetPxPyPzE(lepton1Vector->Px(),lepton1Vector->Py(),lepton1Vector->Pz(),lepton1Vector->E());
      fwlite::Handle< int > l1Flavor;
      l1Flavor.getByLabel(ev,"edmNtuple","lepton1Flavor");
      fwlite::Handle< int > l1Charge;
      l1Charge.getByLabel(ev,"edmNtuple","lepton1Charge");
      fwlite::Handle< ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > lepton2Vector;
      lepton2Vector.getByLabel(ev,"edmNtuple","lepton2LorentzVector");
      TLorentzVector lepton2; lepton2.SetPxPyPzE(lepton2Vector->Px(),lepton2Vector->Py(),lepton2Vector->Pz(),lepton2Vector->E());
      fwlite::Handle<int > l2Flavor;
      l2Flavor.getByLabel(ev,"edmNtuple","lepton2Flavor");
      fwlite::Handle<int > l2Charge;
      l2Charge.getByLabel(ev,"edmNtuple","lepton2Charge");

      
      //cout<<"b-jets"<<endl;


      //b-jets
      fwlite::Handle< ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > bJet1Vector;
      bJet1Vector.getByLabel(ev,"edmNtuple","bJet1LorentzVector");
      TLorentzVector bJet1; bJet1.SetPxPyPzE(bJet1Vector->Px(),bJet1Vector->Py(),bJet1Vector->Pz(),bJet1Vector->E());
      fwlite::Handle< ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > bJet2Vector;
      bJet2Vector.getByLabel(ev,"edmNtuple","bJet2LorentzVector");
      TLorentzVector bJet2; bJet2.SetPxPyPzE(bJet2Vector->Px(),bJet2Vector->Py(),bJet2Vector->Pz(),bJet2Vector->E());
      fwlite::Handle< int > loosebJets;
      loosebJets.getByLabel(ev,"edmNtuple","nLoosebJets");
      fwlite::Handle< int > tightbJets;
      tightbJets.getByLabel(ev,"edmNtuple","nTightbJets");


      //cout<<"MET"<<endl;


      //MET      
      fwlite::Handle< ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > theMET;
      theMET.getByLabel(ev,"edmNtuple","METLorentzVector");
      TLorentzVector METLorentzVector; METLorentzVector.SetPxPyPzE(theMET->Px(),theMET->Py(),theMET->Pz(),theMET->E());


      //cout<<"filling branches"<<endl;


      //fill branches
      lightJet1Pt  = lightJet1.Pt();
      lightJet1Eta = lightJet1.Eta();
      lightJet1Phi = lightJet1.Phi();

      lightJet2Pt  = lightJet2.Pt();
      lightJet2Eta = lightJet2.Eta();
      lightJet2Phi = lightJet2.Phi();

      lepton1Pt     = lepton1.Pt();
      lepton1Eta    = lepton1.Eta();
      lepton1Phi    = lepton1.Phi();
      lepton1Flavor = l1Flavor.ref();
      lepton1Charge = l1Charge.ref();

      lepton2Pt     = lepton2.Pt();
      lepton2Eta    = lepton2.Eta();
      lepton2Phi    = lepton2.Phi();
      lepton2Flavor = l2Flavor.ref();
      lepton2Charge = l2Charge.ref();

      dileptonInvariantMass = (lepton1+lepton2).M();
      sameFlavorLeptons = (abs(lepton1Flavor)==abs(lepton2Flavor));
      inZwindow = (dileptonInvariantMass>76. && dileptonInvariantMass<106.);
      double wideZwindow = (dileptonInvariantMass>66. && dileptonInvariantMass<116.);
      nElectrons=0;
      if (abs(lepton1Flavor)==11) {
	nElectrons++;
	ele1Pt  = lepton1Pt;
	ele1Eta = lepton1Eta;
      }
      if (abs(lepton2Flavor)==11) {
	nElectrons++;
	if (nElectrons==1) {
	  ele1Pt  = lepton2Pt;
	  ele1Eta = lepton2Eta;
	}
      }
      nMuons=0;
      if (abs(lepton1Flavor)==13) {
	nMuons++;
	muon1Pt  = lepton1Pt;
	muon1Eta = lepton1Eta;
      }
      if (abs(lepton2Flavor)==13) {
	nMuons++;
	if (nMuons==1) {
	  muon1Pt  = lepton2Pt;
	  muon1Eta = lepton2Eta;
	}
      }
      doubleEle = (sameFlavorLeptons) ? (nElectrons>1) : 0;
      doubleMu = (sameFlavorLeptons) ? (nMuons>1) : 0;

      bJet1Pt  = bJet1.Pt();
      bJet1Eta = bJet1.Eta();
      bJet1Phi = bJet1.Phi();

      bJet2Pt  = bJet2.Pt();
      bJet2Eta = bJet2.Eta();
      bJet2Phi = bJet2.Phi();

      nLoosebJets = loosebJets.ref();
      nTightbJets = tightbJets.ref();

      MET = METLorentzVector.Et();

      Zveto = (!inZwindow || (inZwindow && !sameFlavorLeptons));
      wideInvertZveto = (wideZwindow && sameFlavorLeptons);
      invertZveto = (inZwindow && sameFlavorLeptons);

      //detector acceptance
      //jet pT cuts chosen to kill QCD
      //invariant mass cut built in for DY sample (M-50)
      kinCut = (abs(lightJet1Eta)<2.4 && abs(lightJet2Eta)<2.4 && abs(bJet1Eta)<2.4 && abs(bJet2Eta)<2.4 && abs(lepton1Eta)<2.4 && abs(lepton2Eta)<2.4 &&
		lightJet2Pt>20 && bJet2Pt>20 && lepton2Pt>10 &&
		dileptonInvariantMass>50);

      tree->Fill();


      if (isRealData) continue; //not looking at data yet

      //top truth information
      fwlite::Handle< bool > truth11;
      truth11.getByLabel(ev,"edmNtuple","top11Truth");
      top11Truth = truth11.ref();
      fwlite::Handle< bool > truth22;
      truth22.getByLabel(ev,"edmNtuple","top22Truth");
      top22Truth = truth22.ref();
      fwlite::Handle< bool > truth12;
      truth12.getByLabel(ev,"edmNtuple","top12Truth");
      top12Truth = truth12.ref();
      fwlite::Handle< bool > truth21;
      truth21.getByLabel(ev,"edmNtuple","top21Truth");
      top21Truth = truth21.ref();

      //sbottom truth information
      fwlite::Handle< bool > truth111;
      truth111.getByLabel(ev,"edmNtuple","sbottomTop11Light1Truth");
      sbottomTop11Light1Truth = truth111.ref();
      fwlite::Handle< bool > truth112;
      truth112.getByLabel(ev,"edmNtuple","sbottomTop11Light2Truth");
      sbottomTop11Light2Truth = truth112.ref();
      fwlite::Handle< bool > truth121;
      truth121.getByLabel(ev,"edmNtuple","sbottomTop12Light1Truth");
      sbottomTop12Light1Truth = truth121.ref();
      fwlite::Handle< bool > truth122;
      truth122.getByLabel(ev,"edmNtuple","sbottomTop12Light2Truth");
      sbottomTop12Light2Truth = truth122.ref();
      fwlite::Handle< bool > truth211;
      truth211.getByLabel(ev,"edmNtuple","sbottomTop21Light1Truth");
      sbottomTop21Light1Truth = truth211.ref();
      fwlite::Handle< bool > truth212;
      truth212.getByLabel(ev,"edmNtuple","sbottomTop21Light2Truth");
      sbottomTop21Light2Truth = truth212.ref();
      fwlite::Handle< bool > truth221;
      truth221.getByLabel(ev,"edmNtuple","sbottomTop22Light1Truth");
      sbottomTop22Light1Truth = truth221.ref();
      fwlite::Handle< bool > truth222;
      truth222.getByLabel(ev,"edmNtuple","sbottomTop22Light2Truth");
      sbottomTop22Light2Truth = truth222.ref();

      //reconstructed tops
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > top11;
      top11.getByLabel(ev,"edmNtuple","top11Vectors");
      top11Vector = top11.ref();
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > top12;
      top12.getByLabel(ev,"edmNtuple","top12Vectors");
      top12Vector = top12.ref();
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > top21;
      top21.getByLabel(ev,"edmNtuple","top21Vectors");
      top21Vector = top21.ref();
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > top22;
      top22.getByLabel(ev,"edmNtuple","top22Vectors");
      top22Vector = top22.ref();


      //reoncstructed sbottoms
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > sbottomTop11Light1;
      sbottomTop11Light1.getByLabel(ev,"edmNtuple","sbottomTop11Light1");
      sbottomTop11Light1Vector = sbottomTop11Light1.ref();
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > sbottomTop11Light2;
      sbottomTop11Light2.getByLabel(ev,"edmNtuple","sbottomTop11Light2");
      sbottomTop11Light2Vector = sbottomTop11Light2.ref();
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > sbottomTop12Light1;
      sbottomTop12Light1.getByLabel(ev,"edmNtuple","sbottomTop12Light1");
      sbottomTop12Light1Vector = sbottomTop12Light1.ref();
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > sbottomTop12Light2;
      sbottomTop12Light2.getByLabel(ev,"edmNtuple","sbottomTop12Light2");
      sbottomTop12Light2Vector = sbottomTop12Light2.ref();
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > sbottomTop21Light1;
      sbottomTop21Light1.getByLabel(ev,"edmNtuple","sbottomTop21Light1");
      sbottomTop21Light1Vector = sbottomTop21Light1.ref();
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > sbottomTop21Light2;
      sbottomTop21Light2.getByLabel(ev,"edmNtuple","sbottomTop21Light2");
      sbottomTop21Light2Vector = sbottomTop21Light2.ref();
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > sbottomTop22Light1;
      sbottomTop22Light1.getByLabel(ev,"edmNtuple","sbottomTop22Light1");
      sbottomTop22Light1Vector = sbottomTop22Light1.ref();
      fwlite::Handle< std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > > sbottomTop22Light2;
      sbottomTop22Light2.getByLabel(ev,"edmNtuple","sbottomTop22Light2");
      sbottomTop22Light2Vector = sbottomTop22Light2.ref();


      //reconstruction selection
      sbottomVectors();

      deltaM = abs(sbottom1.M()-sbottom2.M());

      if (!isRealData) {
      	sbottomCandMass = sbottom1.M();
	sbottomTruth = sbottom1Truth;
      	sbottomTree->Fill();
      	sbottomCandMass = sbottom2.M();
	sbottomTruth = sbottom2Truth;
	sbottomTree->Fill();
      }

      
    } //end loop over events

  }

  else {
    cout<<"I couldn't find the input file"<<endl;
    return;
  }

  fout->Write();
  fout->Close();
  fin->Close();

}

void makeTree::makeTopPairs()
{
  topPairs.clear();
  topPairs.push_back( make_pair(top11Vector,top22Vector) );
  topPairs.push_back( make_pair(top11Vector,top22Vector) );
  topPairs.push_back( make_pair(top12Vector,top21Vector) );
  topPairs.push_back( make_pair(top12Vector,top21Vector) );

  topPairTruths.clear();
  topPairTruths.push_back( make_pair(top11Truth,top22Truth) );
  topPairTruths.push_back( make_pair(top11Truth,top22Truth) );
  topPairTruths.push_back( make_pair(top12Truth,top21Truth) );
  topPairTruths.push_back( make_pair(top12Truth,top21Truth) );
}

void makeTree::makeSbottomPairs()
{
  sbottomPairs.clear();
  sbottomPairs.push_back( make_pair(sbottomTop11Light1Vector,sbottomTop22Light2Vector) );
  sbottomPairs.push_back( make_pair(sbottomTop11Light2Vector,sbottomTop22Light1Vector) );
  sbottomPairs.push_back( make_pair(sbottomTop12Light1Vector,sbottomTop21Light2Vector) );
  sbottomPairs.push_back( make_pair(sbottomTop12Light2Vector,sbottomTop21Light1Vector) );

  sbottomPairTruths.clear();
  sbottomPairTruths.push_back( make_pair(sbottomTop11Light1Truth,sbottomTop22Light2Truth) );
  sbottomPairTruths.push_back( make_pair(sbottomTop11Light2Truth,sbottomTop22Light1Truth) );
  sbottomPairTruths.push_back( make_pair(sbottomTop12Light1Truth,sbottomTop21Light2Truth) );
  sbottomPairTruths.push_back( make_pair(sbottomTop12Light2Truth,sbottomTop21Light1Truth) );
}

void makeTree::makeLightJetPairs()
{
  lightJetPairs.clear();
  lightJetPairs.push_back( make_pair(lightJet1,lightJet2) );
  lightJetPairs.push_back( make_pair(lightJet2,lightJet1) );
  lightJetPairs.push_back( make_pair(lightJet1,lightJet2) );
  lightJetPairs.push_back( make_pair(lightJet2,lightJet1) );
}

void makeTree::initializeStuff()
{
  top1Truth=false; sbottom1Truth=false;
  top2Truth=false; sbottom2Truth=false;
  foundPair=false;
  top1.SetPxPyPzE(0.,0.,0.,-1.); sbottom1.SetPxPyPzE(0.,0.,0.,-1.);
  top2.SetPxPyPzE(0.,0.,0.,-1.); sbottom2.SetPxPyPzE(0.,0.,0.,-1.);
  deltaM=99999.;
  //cout<<"top 1 mass: "<<top1.M()<<endl;
  //cout<<"top 2 mass: "<<top2.M()<<endl;
  //cout<<"sbottom 1 mass: "<<sbottom1.M()<<endl;
  //cout<<"sbottom 2 mass: "<<sbottom2.M()<<endl;
  //cout<<"deltaM: "<<deltaM<<endl;
}

void makeTree::sbottomVectors()
{
  initializeStuff();
  //makeLightJetPairs();
  makeTopPairs();
  makeSbottomPairs();

  double dM = 99999.;
  unsigned int minDeltaMIndex=-1, minDeltaMNeutrinoIndex=-1;
  bool found_dM = false;

  for (unsigned int iPair=0; iPair<sbottomPairs.size(); iPair++) {
    std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > thisSbottom1Vector = (sbottomPairs[iPair]).first;
    std::vector<ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > > thisSbottom2Vector = (sbottomPairs[iPair]).second;
    for (unsigned int iSol=0; iSol<thisSbottom1Vector.size(); iSol++) {
      ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > thisSbottom1 = thisSbottom1Vector[iSol];
      ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > thisSbottom2 = thisSbottom2Vector[iSol];
      if (thisSbottom1.mass()==0 || thisSbottom2.mass()==0) {
        continue;
      }
      double thisDeltaM = abs(thisSbottom1.mass()-thisSbottom2.mass());
      if (thisDeltaM<dM) {
        dM = thisDeltaM;
        minDeltaMIndex = iPair;
        minDeltaMNeutrinoIndex = iSol;
        found_dM = true;
      }
    }
  }
  if (found_dM) {
    ROOT::Math::LorentzVector<ROOT::Math::PxPyPzE4D<double> > v1, v2; 
    TLorentzVector vl1, vl2; 
    bool t1, t2;

    v1 = ((sbottomPairs[minDeltaMIndex]).first )[minDeltaMNeutrinoIndex]; vl1.SetPxPyPzE( v1.Px(), v1.Py(), v1.Pz(), v1.E() );
    v2 = ((sbottomPairs[minDeltaMIndex]).second)[minDeltaMNeutrinoIndex]; vl2.SetPxPyPzE( v2.Px(), v2.Py(), v2.Pz(), v2.E() );
    vl1.SetPxPyPzE( v1.Px(), v1.Py(), v1.Pz(), v1.E() );
    vl2.SetPxPyPzE( v2.Px(), v2.Py(), v2.Pz(), v2.E() );
    sbottom1 = (vl1.M()>vl2.M()) ? vl1 : vl2;
    sbottom2 = (vl1.M()>vl2.M()) ? vl2 : vl1;

    v1 = ((topPairs[minDeltaMIndex]).first )[minDeltaMNeutrinoIndex];
    v2 = ((topPairs[minDeltaMIndex]).second)[minDeltaMNeutrinoIndex];
    vl1.SetPxPyPzE( v1.Px(), v1.Py(), v1.Pz(), v1.E() );
    vl2.SetPxPyPzE( v2.Px(), v2.Py(), v2.Pz(), v2.E() );
    top1 = (sbottom1.M()>sbottom2.M()) ? vl1 : vl2;
    top2 = (sbottom1.M()>sbottom2.M()) ? vl2 : vl1;

    t1 = (sbottomPairTruths[minDeltaMIndex]).first;
    t2 = (sbottomPairTruths[minDeltaMIndex]).second;
    sbottom1Truth = (sbottom1.M()>sbottom2.M()) ? t1 : t2;
    sbottom2Truth = (sbottom1.M()>sbottom2.M()) ? t2 : t1;

    t1 = (topPairTruths[minDeltaMIndex]).first;
    t2 = (topPairTruths[minDeltaMIndex]).second;
    top1Truth = (sbottom1.M()>sbottom2.M()) ? t1 : t2;
    top2Truth = (sbottom1.M()>sbottom2.M()) ? t2 : t1;

    foundPair=true;
  }
  else {
    TLorentzVector zero(0,0,0,-1);
    sbottom1=zero;
    sbottom2=zero;
    top1=zero;
    top2=zero;
    sbottom1Truth=false;
    sbottom2Truth=false;
    foundPair=false;
  }
  //cout<<"top 1 mass: "<<top1.M()<<endl;
  //cout<<"top 2 mass: "<<top2.M()<<endl;
  //if (foundPair) {
  //  cout<<"sbottom 1 mass: "<<sbottom1.M()<<", truth: "<<boolalpha<<sbottom1Truth<<endl;
  //  cout<<"sbottom 2 mass: "<<sbottom2.M()<<", truth: "<<boolalpha<<sbottom2Truth<<endl;
  //  //cout<<"deltaM: "<<dM<<endl;
  //}
  return;
}