//gROOT->Reset();
#include "strCut.h"
#include "TMath.h"
#include "TBranch.h"
#include "TTree.h"
#include "TFile.h"
#include "TH1F.h"
#include <fstream>
#include "TH2F.h"
#include <sstream>
#include "string.h"
#include <iostream>
#include <cmath>

double highest_pt[6][9];
enum ENUM_List{highpt_pt,highpt_px,highpt_py,highpt_pz,highpt_e,highpt_eta,highpt_phi,highpt_charge,highpt_index};
enum ENUM_List2{allelec,looseelec,medelec,tightelec,allmu,combmu};
enum ENUM_List3{m_all_e,m_loose_e,m_med_e,m_tight_e};
enum ENUM_List4{m_all_mu,m_comb_mu};

int n_each_type[6]={0,0,0,0,0,0};

double phiDif(double phi1,double phi2){
  double pi = TMath::Pi();
  double dif = fabs(phi1-phi2);
  if (dif > pi) {
    if (phi1>phi2){dif = phi1 - (phi2 + 2*pi);} else {dif = phi2 - (phi1 + 2*pi);}
    dif = fabs(dif);
  }
  return dif;
}

double cInvMass(double eE,double ePx,double ePy,double ePz,double mE,double mPx,double mPy,double mPz ){
  double mP_dot_eP = (ePx*mPx + ePy*mPy + ePz*mPz);
  return sqrt( eE*eE + mE*mE - ePx*ePx - ePy*ePy - ePz*ePz - mPx*mPx - mPy*mPy - mPz*mPz + 2*( eE*mE - mP_dot_eP) );
}

void highest_pt_extract_mass ( int particle, double charge, double pt, double px, double py, double pz, double e, double phi, double eta, int index){
	if (pt > highest_pt[particle][highpt_pt]){
		highest_pt[particle][highpt_pt]=pt;
		highest_pt[particle][highpt_px]=px;
		highest_pt[particle][highpt_py]=py;
		highest_pt[particle][highpt_pz]=pz;
		highest_pt[particle][highpt_e]=e;
		highest_pt[particle][highpt_eta]=eta;
		highest_pt[particle][highpt_phi]=phi;
		highest_pt[particle][highpt_charge]=charge;
		highest_pt[particle][highpt_index]=index;
	}
	return;
}


void readEmu(string strInFilePrefix,int nFiles, string strOutFilePrefix = "analysis",string strBranchNameTxt = "branchNames_private.txt"){
 
  //TFile *f =new TFile(strInFile.c_str(), "READ");
  //TTree *tree = (TTree*)gDirectory->Get("physics");

  Int_t         el_n;
  vector<float> *el_e;
  vector<float> *el_pt;
  vector<float> *el_px;
  vector<float> *el_py;
  vector<float> *el_pz;
  vector<float> *el_eta;
  vector<float> *el_phi;
  vector<float> *el_charge;
  vector<int>   *el_loose;
  vector<int>   *el_medium;
  vector<int>   *el_tight;

  Int_t         mu_n;
  vector<float> *mu_e;
  vector<float> *mu_pt;
  vector<float> *mu_px;
  vector<float> *mu_py;
  vector<float> *mu_pz;
  vector<float> *mu_eta;
  vector<float> *mu_phi;
  vector<float> *mu_charge;
  vector<int>   *mu_isCombinedMuon;
  vector<float> *mu_etcone40;
  vector<int>   *mu_bestMatch;

  float MET_et;
  Int_t jet_n;
  vector<float> *jet_pt;
  Bool_t         EF_mu4;
  Bool_t         EF_e10_medium;
  Bool_t         EF_mu6;
  Bool_t         EF_mu10;

  // List of branches
  TBranch        *b_el_n; 
  TBranch        *b_el_e;
  TBranch        *b_el_pt;
  TBranch        *b_el_px;
  TBranch        *b_el_py;
  TBranch        *b_el_pz;
  TBranch        *b_el_eta;
  TBranch        *b_el_phi;
  TBranch        *b_el_charge;   
  TBranch        *b_el_loose;   
  TBranch        *b_el_medium;   
  TBranch        *b_el_tight;   

  TBranch        *b_mu_n;  
  TBranch        *b_mu_e; 
  TBranch        *b_mu_pt; 
  TBranch        *b_mu_px; 
  TBranch        *b_mu_py; 
  TBranch        *b_mu_pz; 
  TBranch        *b_mu_eta;  
  TBranch        *b_mu_phi;  
  TBranch        *b_mu_charge;   
  TBranch        *b_mu_isCombinedMuon;   
  TBranch        *b_mu_etcone40;
  TBranch        *b_mu_bestMatch;

  TBranch        *b_MET_et;  
  TBranch        *b_jet_n;  
  TBranch        *b_jet_pt; 

  TBranch         *b_EF_mu4;
  TBranch         *b_EF_e10_medium;
  TBranch         *b_EF_mu6;
  TBranch         *b_EF_mu10;
  
  TH1F *h_invmass_comb[combo_n_Hist];
  TH1F *h_invmass_comb_opp_charge[combo_n_Hist];
  TH1F *h_invmass_comb_same_charge[combo_n_Hist];
  TH1F *h_MET_comb[combo_n_Hist];
  TH1F *h_phi_diff_comb[combo_n_Hist];
  TH1F *h_invmass_comb_pt_cut[combo_n_Hist];
  TH1F *h_invmass_comb_MET_cut[combo_n_Hist];
  TH1F *h_invmass_comb_dphi_cut[combo_n_Hist];
  TH1F *h_invmass_comb_all_cut[combo_n_Hist];
  TH1F *h_lepton_pt_comb[combo_n_Hist];
  TH2F *h_invmass_v_MET_comb[combo_n_Hist];
  TH2F *h_invmass_v_dphi_comb[combo_n_Hist];
  TH1F *h_invmass_low_dphi_comb[combo_n_Hist];
  TH1F *h_cuts;
  TH1F *h_cutsPercent;
  TFile *f;
  TTree *tree;
  
  
  TH1F *h_jetn = new TH1F("h_jetn","Number of Jets after other cuts",26,-0.5,25.5);
  TH1F *h_phiDif = new TH1F("h_phiDif","Phi Difference between e #mu pair",100,0,3.2);
  /*TH1F *h_muPtOrder = new TH1F("h_muPtOrder","Highest Pt Muon Index",totHist+1,-0.5,(double)totHist+0.5);
  TH1F *h_muPtOrder_combined = new TH1F("h_muPtOrder_combined","Highest Pt Muon Index",totHist+1,-0.5,(double)totHist+0.5);
  TH1F *h_elPtOrder = new TH1F("h_elPtOrder","Highest Pt Electron Index",totHist+1,-0.5,(double)totHist+0.5);
  TH1F *h_elPtOrder_loose = new TH1F("h_elPtOrder_loose","Highest Pt Electron Index",totHist+1,-0.5,(double)totHist+0.5);
  TH1F *h_elPtOrder_med = new TH1F("h_elPtOrder_med","Highest Pt Electron Index",totHist+1,-0.5,(double)totHist+0.5);
  TH1F *h_elPtOrder_tight = new TH1F("h_elPtOrder_tight","Highest Pt Electron Index",totHist+1,-0.5,(double)totHist+0.5);*/
  TH1F *h_muonCut = new TH1F("h_muonCut","Muon Cut",3,-0.5,2.5);
  TH1F *h_MET = new TH1F("h_MET","Missing ET",100,0,1000000);
  
  TH1F *h_trigEff;
  
  char buffTitle[200],buffCut[200];

  /*tree->SetBranchAddress("el_n", &el_n, &b_el_n);
  tree->SetBranchAddress("el_E", &el_e, &b_el_e);
  tree->SetBranchAddress("el_pt", &el_pt, &b_el_pt);
  tree->SetBranchAddress("el_px", &el_px, &b_el_px);
  tree->SetBranchAddress("el_py", &el_py, &b_el_py);
  tree->SetBranchAddress("el_pz", &el_pz, &b_el_pz);
  tree->SetBranchAddress("el_eta", &el_eta, &b_el_eta);
  tree->SetBranchAddress("el_phi", &el_phi, &b_el_phi);
  tree->SetBranchAddress("el_charge", &el_charge, &b_el_charge);
  tree->SetBranchAddress("el_loose", &el_loose, &b_el_loose);
  tree->SetBranchAddress("el_medium", &el_medium, &b_el_medium);
  tree->SetBranchAddress("el_tight", &el_tight, &b_el_tight);

  tree->SetBranchAddress("mu_n", &mu_n, &b_el_n);
  tree->SetBranchAddress("mu_E", &mu_e, &b_mu_e);
  tree->SetBranchAddress("mu_pt", &mu_pt, &b_mu_pt);
  tree->SetBranchAddress("mu_px", &mu_px, &b_mu_px);
  tree->SetBranchAddress("mu_py", &mu_py, &b_mu_py);
  tree->SetBranchAddress("mu_pz", &mu_pz, &b_mu_pz);
  tree->SetBranchAddress("mu_eta", &mu_eta, &b_mu_eta);
  tree->SetBranchAddress("mu_phi", &mu_phi, &b_mu_phi);
  tree->SetBranchAddress("mu_charge", &mu_charge, &b_mu_charge);
  tree->SetBranchAddress("mu_isCombinedMuon", &mu_isCombinedMuon, &b_mu_isCombinedMuon);
  tree->SetBranchAddress("mu_etcone40", &mu_etcone40, &b_mu_etcone40);
  tree->SetBranchAddress("mu_bestMatch", &mu_bestMatch, &b_mu_bestMatch);

  tree->SetBranchAddress("MET_et", &MET_et, &b_MET_et);
  tree->SetBranchAddress("jet_n", &jet_n, &b_jet_n);

  tree->SetBranchAddress("EF_mu4", &EF_mu4, &b_EF_mu4);
  tree->SetBranchAddress("EF_e10_medium", &EF_e10_medium, &b_EF_e10_medium);
  tree->SetBranchAddress("EF_mu6", &EF_mu6, &b_EF_mu6);
  tree->SetBranchAddress("EF_mu10", &EF_mu10, &b_EF_mu10);*/
  
  for (int iFiles = 0; iFiles < nFiles; iFiles++){ 
    char buffInFile[400];
    if (nFiles ==1) sprintf(buffInFile,"%s",strInFilePrefix.c_str());
    else if (nFiles >= 10) {
      sprintf(buffInFile,"%s%.2i.root",strInFilePrefix.c_str(),iFiles+1);
    }
    else {
     sprintf(buffInFile,"%s%i.root",strInFilePrefix.c_str(),iFiles+1);
    }
    

    f =new TFile(buffInFile, "READ");
    tree = (TTree*)gDirectory->Get("physics");
    ifstream bNames (strBranchNameTxt.c_str());
    string strBNames[50];
    int bNameCount=0;
    tree->SetBranchStatus("*",0);
    while (!bNames.eof()){
      getline (bNames,strBNames[bNameCount]);
      if (strBNames[bNameCount]!=""){
      tree->SetBranchStatus(strBNames[bNameCount].c_str(),1);}
      bNameCount++;
    }
    bNameCount=0;
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_n, &b_el_n);                    
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_e, &b_el_e);                    
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_pt, &b_el_pt);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_px, &b_el_px);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_py, &b_el_py);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_pz, &b_el_pz);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_eta, &b_el_eta);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_phi, &b_el_phi);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_charge, &b_el_charge);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_loose, &b_el_loose);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_medium, &b_el_medium);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &el_tight, &b_el_tight);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_n, &b_mu_n);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_e, &b_mu_e);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_pt, &b_mu_pt);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_px, &b_mu_px);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_py, &b_mu_py);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_pz, &b_mu_pz);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_eta, &b_mu_eta);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_phi, &b_mu_phi);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_charge, &b_mu_charge);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_isCombinedMuon, &b_mu_isCombinedMuon);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_etcone40, &b_mu_etcone40);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &mu_bestMatch, &b_mu_bestMatch);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &MET_et, &b_MET_et);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &jet_n, &b_jet_n);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &jet_pt, &b_jet_pt);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &EF_mu4, &b_EF_mu4);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &EF_e10_medium, &b_EF_e10_medium);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &EF_mu6, &b_EF_mu6);
    tree->SetBranchAddress(strBNames[bNameCount++].c_str(), &EF_mu10, &b_EF_mu10);


  Long64_t nentries = tree->GetEntriesFast();

  cout << "Total events: " << nentries << endl;
  char buffTitle[200],buffCut[200];

  TH1F *h_invmass_cut[totHist];
  string strCTitle = "Events vs cuts;0=Events ";
  for (int j = 0; j < totHist; j++){
    //It looks messy, but c++ doesn't do strings well and it works
    h_invmass_cut[j] = new TH1F(((string)("h_invmass_cut_"+strCut[j])).c_str(),((string)("e #mu Invariant Mass "+strCutTitle[j]+"; Invariant Mass (MeV)")).c_str(),binN,0,maxBin);
    std::ostringstream oss; oss << (j+1); //The lack of ability to handle string conversions in C++ is infuriating
    strCTitle += oss.str() + "=" + strCut[j] + " "; 
  }
  
  //string strCombTitle = "Events vs combos; ";
  for (int j = 0; j < combo_n_Hist; j++){
    //It looks messy, but c++ doesn't do strings well and it works
    h_invmass_comb[j] = new TH1F(((string)("h_invmass_comb_"+strComb[j])).c_str(),((string)("e #mu Invariant Mass "+strCombTitle[j]+"; Invariant Mass (MeV)")).c_str(),binN,0,maxBin);
    //string testy=("h_invmass_comb_"+strComb[j]);
    //string testy2=("e #mu Invariant Mass "+strCombTitle[j]+"; Invariant Mass (MeV)");
    //string testy3=(testy2+"; Invariant Mass (MeV)");
    /*char tmp_buff[300];
    char tmp_buff2[300];
    sprintf(tmp_buff2,"e mu Invariant Mass %s Invariant Mass (MeV)",strComb[j].c_str());
    sprintf(tmp_buff,"h_invmass_comb_%s",strCombTitle[j].c_str());
    h_invmass_comb[j] = new TH1F(tmp_buff,tmp_buff2,binN,0,maxBin);*/
    h_invmass_comb_opp_charge[j] = new TH1F(((string)("h_invmass_comb_opp_charge"+strComb[j])).c_str(),((string)("e #mu Invariant Mass (opp charge)"+strCombTitle[j]+"; Invariant Mass (MeV)")).c_str(),binN,0,maxBin);
    h_invmass_comb_same_charge[j] = new TH1F(((string)("h_invmass_comb_same_charge"+strComb[j])).c_str(),((string)("e #mu Invariant Mass (same charge)"+strCombTitle[j]+"; Invariant Mass (MeV)")).c_str(),binN,0,maxBin);
    h_invmass_comb_all_cut[j] = new TH1F(((string)("h_invmass_comb_cut_"+strComb[j])).c_str(),((string)("e #mu Invariant Mass "+strCombTitle[j]+" with Pt, MET, opp charge and Phi cuts; Invariant Mass (MeV)")).c_str(),binN,0,maxBin);
    h_MET_comb[j] = new TH1F(((string)("h_MET_comb_"+strComb[j])).c_str(),((string)("e #mu event Missing ET "+strCombTitle[j]+"; MET (MeV)")).c_str(),binN,0,maxBin);
    h_lepton_pt_comb[j] = new TH1F(((string)("h_lepton_pt_comb_"+strComb[j])).c_str(),((string)("e #mu Lepton p_{T} "+strCombTitle[j]+"; p_{T} (MeV)")).c_str(),binN,0,maxBin);
    h_phi_diff_comb[j] = new TH1F(((string)("h_phi_diff_comb_"+strComb[j])).c_str(),((string)("e #mu #phi difference "+strCombTitle[j]+"; Phi (Radians)")).c_str(),100,0,3.2);
    h_invmass_v_MET_comb[j] = new TH2F(((string)("h_invmass_v_MET_comb_"+strComb[j])).c_str(),((string)("e #mu Inv Mass vs. MET "+strCombTitle[j]+"; InvMass (MeV); MET (MeV)")).c_str(),50,0,2000000,50,0,2000000);
    h_invmass_v_dphi_comb[j] = new TH2F(((string)("h_invmass_v_dphi_comb_"+strComb[j])).c_str(),((string)("e #mu Inv Mass vs. dPhi "+strCombTitle[j]+"; Inv Mass(MeV); Phi (Radians)")).c_str(),50,0,2000000,50,0,3.2);
    h_invmass_low_dphi_comb[j] = new TH1F(((string)("h_invmass_low_dphi_comb_"+strComb[j])).c_str(),((string)("e #mu Inv Mass at Low #phi difference "+strCombTitle[j]+"; (MeV)")).c_str(),100,0,2000000);
    //std::ostringstream oss; oss << (j+1); //The lack of ability to handle string conversions in C++ is infuriating
    //strCombTitle = strCombTitle + oss.str() + "=" + strComb[j] + " "; 
  }
  cout << strCTitle << endl;

  h_cuts = new TH1F("h_cuts",strCTitle.c_str(),totHist+1,-0.5,(double)totHist+0.5);
  h_cutsPercent = new TH1F("h_cutsPercent",strCTitle.c_str(),11,-0.5,10.5);

  int trigCount[trigN];
  for (int i = 0; i < trigN; i++) trigCount[i] = 0;
  string strTTitle = "Trigger Efficiencies;";
  for (int j = 0; j < trigN; j++){
    //It looks messy, but c++ doesn't do strings well and it works
    std::ostringstream oss; oss << (j); //The lack of ability to handle string conversions in C++ is infuriating
    strTTitle += oss.str() + "=" + strTrigs[j] + " "; 
  }
  cout << strTTitle << endl;
  
  TH1F *h_trigEff = new TH1F("h_trigEff",strTTitle.c_str(),12,-0.5,11.5);

  int count = 0;
  int trigEfCount=0;
  for (int i = 0; i < nentries; i++){
    tree->GetEntry(i);
   
	
    if(count++ == nentries / 10 ) {cout << "Event #: " << i << " (" << abs(100*((double)i/(double)nentries)) << "% processed)" << endl; count=0;}  //Prints out the current event number every ~1/10 of the way through running
    if (el_n > 0 && mu_n > 0 ){ 
      /*double mu_highPt[2] = {0,0}; //Index: 0 = no quality cut, 1 = combined muons
      double el_highPt[4]= {0,0,0,0}; //Index: 0 = no quality cut, 1 = loose, 2 = medium, 3 = tight
      double ePt[4], eCharge[4], ePx[4], ePy[4], ePz[4], eE[4], ePhi[4], eEta[4], eHighPtIndex[4]; //Index: 0 = no quality cut, 1 = loose, 2 = medium, 3 = tight
      double mPt[2], mCharge[2], mPx[2], mPy[2], mPz[2], mE[2], mPhi[2], mEta[2], mHighPtIndex[2],mBestMatch[2],mEtCone40[2]; //Index: 0 = no quality cut, 1 = combined muons*/
      
      for (int k=0; k<6; k++){
        n_each_type[k]=0;
      	for (int j=0; j<9;j++){
		highest_pt[k][j]=0;
	}
      }


      //********************************
      // Get highest Pt all/loose/medium/tight electrons
      //********************************

      for (int iElect = 0; iElect < el_n; iElect++){
      	if ( fabs(el_eta->at(iElect)) <2.5 && (fabs(el_eta->at(iElect))<1.37 || fabs(el_eta->at(iElect))>1.52)){//elec in ecal
	//All electrons
		highest_pt_extract_mass(allelec,el_charge->at(iElect),el_pt->at(iElect),el_px->at(iElect),el_py->at(iElect),
			el_pz->at(iElect),el_e->at(iElect),el_phi->at(iElect),el_eta->at(iElect),iElect);
		n_each_type[0]++;
		if (el_loose->at(iElect)==1){//loose elestrons
			n_each_type[1]++;
			highest_pt_extract_mass(looseelec,el_charge->at(iElect),el_pt->at(iElect),el_px->at(iElect),
			el_py->at(iElect),el_pz->at(iElect),el_e->at(iElect),el_phi->at(iElect),el_eta->at(iElect),iElect);
		}
		if (el_medium->at(iElect)==1 && el_loose->at(iElect)==1){//med electrons
			highest_pt_extract_mass(medelec,el_charge->at(iElect),el_pt->at(iElect),el_px->at(iElect),
			el_py->at(iElect),el_pz->at(iElect),el_e->at(iElect),el_phi->at(iElect),el_eta->at(iElect),iElect);
			n_each_type[2]++;
		}
		if (el_tight->at(iElect)==1){//tight electrons
			highest_pt_extract_mass(tightelec,el_charge->at(iElect),el_pt->at(iElect),el_px->at(iElect),
			el_py->at(iElect),el_pz->at(iElect),el_e->at(iElect),el_phi->at(iElect),el_eta->at(iElect),iElect);
			n_each_type[3]++;
		}
	}	    
      }
    
      //********************************
      // Get highest Pt all/combined muons
      //********************************
      for (int iMuon = 0; iMuon < mu_n; iMuon++){
      	if ( fabs(mu_eta->at(iMuon))<2.7){//Muon in MS
		highest_pt_extract_mass(allmu,mu_charge->at(iMuon),mu_pt->at(iMuon),mu_px->at(iMuon),mu_py->at(iMuon),
			mu_pz->at(iMuon),mu_e->at(iMuon),mu_phi->at(iMuon),mu_eta->at(iMuon),iMuon);
		n_each_type[4]++;
		if (mu_isCombinedMuon->at(iMuon)==1){//Combined Muons	
			highest_pt_extract_mass(combmu,mu_charge->at(iMuon),mu_pt->at(iMuon),mu_px->at(iMuon),mu_py->at(iMuon),
			mu_pz->at(iMuon),mu_e->at(iMuon),mu_phi->at(iMuon),mu_eta->at(iMuon),iMuon);
			n_each_type[5]++;
		}
	}
      }
      
      /*h_muPtOrder->Fill(highest_pt[allmu][highpt_index]);
      h_muPtOrder_combined->Fill(highest_pt[combmu][highpt_index]);
      h_elPtOrder->Fill(highest_pt[allelec][highpt_index]);
      h_elPtOrder_loose->Fill(highest_pt[looseelec][highpt_index]);
      h_elPtOrder_med->Fill(highest_pt[medelec][highpt_index]);
      h_elPtOrder_tight->Fill(highest_pt[tightelec][highpt_index]);*/

  

      //********************************    
      // Calculate every possible invMass combination with different quality electrons and muons
      //********************************
      
      double invMass[4][2];
      for(int eQ = 0; eQ < 4; eQ++){
	/*invMass[eQ][m_all_mu] = cInvMass( highest_pt[eQ][highpt_e], highest_pt[eQ][highpt_px], highest_pt[eQ][highpt_py],
		highest_pt[eQ][highpt_pz],highest_pt[m_all_mu][highpt_e],highest_pt[m_all_mu][highpt_px],
		highest_pt[m_all_mu][highpt_py], highest_pt[m_all_mu][highpt_pz]);
	invMass[eQ][m_comb_mu] = cInvMass( highest_pt[eQ][highpt_e], highest_pt[eQ][highpt_px], highest_pt[eQ][highpt_py],
		highest_pt[eQ][highpt_pz],highest_pt[m_comb_mu][highpt_e],highest_pt[m_comb_mu][highpt_px],
		highest_pt[m_comb_mu][highpt_py], highest_pt[m_comb_mu][highpt_pz]);*/
	
	int e_index=highest_pt[eQ][highpt_index];
	int mu_index=highest_pt[allmu][highpt_index];
	invMass[eQ][0]=cInvMass(el_e->at(e_index),el_px->at(e_index),el_py->at(e_index),el_pz->at(e_index),mu_e->at(mu_index),
		mu_px->at(mu_index),mu_py->at(mu_index),mu_pz->at(mu_index));
	mu_index=highest_pt[combmu][highpt_index];
	invMass[eQ][1]=cInvMass(el_e->at(e_index),el_px->at(e_index),el_py->at(e_index),el_pz->at(e_index),mu_e->at(mu_index),
		mu_px->at(mu_index),mu_py->at(mu_index),mu_pz->at(mu_index));
	}

       //********************************    
      // Here's where you fill all the combination histograms, without cuts
      // The first index of invMass is electron quality (0=all, 1=loose, 2=medium, 3=tight)
      // The second index is isMuCombined (0=false or true, 1=true)
      // So invMass[1][1] would be the inv mass of the highest pt loose electron and highest pt combined muon
      // The same labeling is used for e & mu variables (ie ePt[1] is the Pt of the highest pt loose electron)
      //********************************
      
      double ptCut = 50*GeV;
      double metCut = 100*GeV;
      for (int w=0; w<4; w++){
      	int e_index=highest_pt[w][highpt_index];
	int mu_index=highest_pt[allmu][highpt_index];
      	if (n_each_type[w]!=0&&n_each_type[4]!=0){
		h_invmass_comb[w]->Fill(invMass[w][0]);
		if ( (el_charge->at(e_index) + mu_charge->at(mu_index))==0){
			h_invmass_comb_opp_charge[w]->Fill(invMass[w][0]);
			h_invmass_v_dphi_comb[w]->Fill(invMass[w][0],phiDif(el_phi->at(e_index),mu_phi->at(mu_index)));
			h_invmass_v_MET_comb[w]->Fill(invMass[w][0],MET_et);
			}
		else{
			h_invmass_comb_same_charge[w]->Fill(invMass[w][0]);
			}
		h_MET_comb[w]->Fill(MET_et);
		h_phi_diff_comb[w]->Fill(phiDif(el_phi->at(e_index),mu_phi->at(mu_index)));
		if ( phiDif(el_phi->at(e_index),mu_phi->at(mu_index))<1){
			h_invmass_low_dphi_comb[w]->Fill(invMass[w][0]);
			}
		h_lepton_pt_comb[w]->Fill(el_pt->at(e_index));
		h_lepton_pt_comb[w]->Fill(mu_pt->at(mu_index));
		}//All Muons
	if (n_each_type[w]!=0&&n_each_type[4]!=0 && MET_et<metCut && el_pt->at(e_index)>ptCut && 
		mu_pt->at(mu_index)>ptCut && phiDif(el_phi->at(e_index),mu_phi->at(mu_index))>2 && 
		(el_charge->at(e_index) + mu_charge->at(mu_index))==0)
		{
		h_invmass_comb_all_cut[w]->Fill(invMass[w][0]);
		}//All Muons with cuts
	if (n_each_type[w]!=0&&n_each_type[5]!=0){
		mu_index=highest_pt[combmu][highpt_index];
		h_invmass_comb[w+4]->Fill(invMass[w][1]);
		if ( (el_charge->at(e_index) + mu_charge->at(mu_index))==0){
			h_invmass_comb_opp_charge[w+4]->Fill(invMass[w][1]);
			h_invmass_v_dphi_comb[w+4]->Fill(invMass[w][1],phiDif(el_phi->at(e_index),mu_phi->at(mu_index)));
			h_invmass_v_MET_comb[w+4]->Fill(invMass[w][1],MET_et);
			}
		else{
			h_invmass_comb_same_charge[w+4]->Fill(invMass[w][1]);
			}
		h_MET_comb[w+4]->Fill(MET_et);
		h_phi_diff_comb[w+4]->Fill(phiDif(el_phi->at(e_index),mu_phi->at(mu_index)));
		if ( phiDif(el_phi->at(e_index),mu_phi->at(mu_index))<1){
			h_invmass_low_dphi_comb[w+4]->Fill(invMass[w][1]);
			}
		h_lepton_pt_comb[w+4]->Fill(el_pt->at(e_index));
		h_lepton_pt_comb[w+4]->Fill(mu_pt->at(mu_index));
		}//Combined Muons
	if (n_each_type[w]!=0&&n_each_type[5]!=0 && MET_et<metCut && el_pt->at(e_index)>ptCut && 
		mu_pt->at(mu_index)>ptCut && phiDif(el_phi->at(e_index),mu_phi->at(mu_index))>2&& 
		(el_charge->at(e_index) + mu_charge->at(mu_index))==0){
		h_invmass_comb_all_cut[w+4]->Fill(invMass[w][0]);
		}//Combined Muons with cuts
	}
      
      
      //********************************    
      // Here's where you fill all the cut histograms
      // The first index of invMass is electron quality (0=all, 1=loose, 2=medium, 3=tight)
      // The second index is isMuCombined (0=false, 1=true)
      // So invMass[1][1] would be the inv mass of the highest pt loose electron and highest pt combined muon
      // The same labeling is used for e & mu variables (ie ePt[1] is the Pt of the highest pt loose electron)
      //********************************
      
      h_invmass_cut[0]->Fill(invMass[0][0]);
      /*if (ePt[0] > ptCut && mPt[0] > ptCut){//pt cut
	h_invmass_cut[1]->Fill(invMass[0][0]);
	
	if (MET_et < metCut){ //MET cut
	  h_invmass_cut[2]->Fill(invMass[0][0]);
	  if (mPt[1] > ptCut && mEtCone40[1] < 1000) {//combined muon cut & isolation cut
	    h_invmass_cut[3]->Fill(invMass[0][1]);
	    if (ePt[1] > ptCut){//loose electron cut
	      h_invmass_cut[4]->Fill(invMass[1][1]);
	      if (ePt[2] > ptCut){//medium electron cut
		h_invmass_cut[5]->Fill(invMass[2][1]);
		if (ePt[3] > ptCut){//tight electron cut
		  h_invmass_cut[6]->Fill(invMass[3][1]); 
		  if (phiDif(ePhi[3],mPhi[1])>2.0){//delta phi > 2.0
		    h_invmass_cut[7]->Fill(invMass[3][1]);
		    if (fabs(eEta[3]) <= 2.5 && (fabs(eEta[3]) >= 1.52 || fabs(eEta[3]) <= 1.37) && fabs(mEta[1]) < 2.7){// muon & electron eta in detector
		      h_invmass_cut[8]->Fill(invMass[3][1]);
		      h_jetn->Fill(jet_n);
		      h_MET->Fill(MET_et);
		      if (jet_n <= 6){
			h_invmass_cut[9]->Fill(invMass[3][1]);
			
			trigEfCount++; //Counts the total number of events to use when calculating trigger efficiencies
			trigCount[0] += EF_mu4;
			trigCount[1] += EF_e10_medium;
			trigCount[2] += EF_mu6;
			if (EF_mu4 && EF_e10_medium) trigCount[3]++;
			if (EF_mu4 || EF_e10_medium) trigCount[4]++;
			if (EF_mu6 && EF_e10_medium) trigCount[5]++;
			if (EF_mu6 || EF_e10_medium) trigCount[6]++;
			trigCount[7] += EF_mu10;
			if (EF_mu10 || EF_e10_medium) {
			  trigCount[8]++;
			  h_invmass_cut[10]->Fill(invMass[3][1]);
			  if (eCharge[3]+mCharge[1] == 0) {
			    h_invmass_cut[11]->Fill(invMass[3][1]);
			  }
			}
			if (EF_mu10 && EF_e10_medium) {
			  trigCount[9]++;
			  //h_invmass_cut[10]->Fill(invMass[3][1]);
			}
		      }
		    }
		  }	    
		} 
	      }	    
	    }
	  }
	}
      }*/
    }
    }
  }

  //********************************    
  // Saves all histograms created from "File.root" to "File_histos.root"
  // Also creates a cut flow histogram (h_cuts) which tracks the total entries in each invMass histogram
  //********************************    



  //int pos = strInFile.find(".root");
  //string strOutFile = strInFile.substr(0,pos) + "_histos.root";
  string strOutFile = strOutFilePrefix + "_histos.root";
  TFile *fout = new TFile(strOutFile.c_str(),"RECREATE");
  //int totEntries = nentries;
  Long64_t nentries = tree->GetEntriesFast();
  h_cuts->SetBinContent(1,nentries);
  h_cutsPercent->SetBinContent(1,1);
  //cout << "Seg fault here?"<<endl;
  /*for (int i = 0; i < totHist; i++){
    int hEntries = h_invmass_cut[i]->GetEntries();
    h_cuts->SetBinContent(i+2,hEntries);
    //cout << i << " | " << hEntries << " | " << h_cuts->GetBinContent(i+1) << endl;
    h_cutsPercent->SetBinContent(i+2,1-((double) hEntries / (double)h_cuts->GetBinContent(i+1)));
    totEntries += hEntries;
    h_invmass_cut[i]->Write();
  }
  h_cuts->SetEntries(totEntries);*/
  
  for (int j=0; j<combo_n_Hist; j++){
  	h_invmass_comb[j]->Write();
	h_MET_comb[j]->Write();
	h_phi_diff_comb[j]->Write();
	h_lepton_pt_comb[j]->Write();
	h_invmass_comb_all_cut[j]->Write();
	h_invmass_comb_opp_charge[j]->Write();
	h_invmass_comb_same_charge[j]->Write();
	h_invmass_v_dphi_comb[j]->Write();
	h_invmass_v_MET_comb[j]->Write();
	h_invmass_low_dphi_comb[j]->Write();
	}

  //cout << "Seg fault here instead?"<<endl;

  cout <<"--------------------" << endl << "Trigger Efficiencies" << endl << "--------------------" << endl;
  /*for (int j = 0; j < trigN; j++){
    h_trigEff->SetBinContent(j+1, (double)trigCount[j] / (double)trigEfCount);//h_cuts->GetBinContent(totHist+1));
    cout << strTrigs[j] << ": " << (double)trigCount[j] / (double)trigEfCount << endl;;//h_cuts->GetBinContent(totHist+1) << endl;
  }*/

  h_cutsPercent->Write();
  h_cuts->Write();
  h_jetn->Write();
  //h_trigEff->Write();

  /*h_elPtOrder->Write();
  h_elPtOrder_loose->Write();
  h_elPtOrder_med->Write();
  h_elPtOrder_tight->Write();
  h_muPtOrder->Write();
  h_muPtOrder_combined->Write();*/
  h_muonCut->Write();
  h_MET->Write();

  f->Close();
  fout->Close();
  cout << "Histograms saved to: " << strOutFile << endl;
}