#include <iostream>
#include <vector>
#include <set>
#include <cmath>
#include <string>
#include <cstdio>
#include <cstring>
#include "TString.h"
#include "TLorentzVector.h"
#include "TVector2.h"
#include "TFile.h"
#include "TTree.h"
#include "TH1.h"
#include "THStack.h"
#include <stdlib.h>
#include "link.h"
using namespace std;

 

int main() {
 TFile f("stopFlatNtuples_1.root");
 TTree *t=(TTree*)f.Get("stopTreeMaker/AUX");


vector<TLorentzVector> *genDecayLVec =0, *jetsLVec=0, *muonsLVec=0;
vector<int> *genDecayIdxVec =0, *genDecayPdgIdVec =0, *genDecayMomIdxVec =0;
vector<string> *genDecayStrVec =0;
vector<int> *W_emuVec=0, *W_tau_emuVec=0, *W_tau_prongsVec=0;
vector<double> *muonsRelIso=0;
int nMuons_CUT=0, nElectrons_CUT=0, nJets_CUT=0;
double ht=0, mht=0, met=0, metphi=0, dPhi0_CUT=0,dPhi1_CUT=0, dPhi2_CUT=0;

t->SetBranchAddress("genDecayIdxVec", &genDecayIdxVec);
t->SetBranchAddress("genDecayPdgIdVec", &genDecayPdgIdVec);
t->SetBranchAddress("genDecayMomIdxVec", &genDecayMomIdxVec);
t->SetBranchAddress("genDecayLVec", &genDecayLVec);
t->SetBranchAddress("genDecayStrVec", &genDecayStrVec);
t->SetBranchAddress("W_emuVec", &W_emuVec);
t->SetBranchAddress("W_tau_emuVec", &W_tau_emuVec);
t->SetBranchAddress("W_tau_prongsVec", &W_tau_prongsVec);
t->SetBranchAddress("jetsLVec", &jetsLVec);
t->SetBranchAddress("muonsLVec", &muonsLVec);
t->SetBranchAddress("muonsRelIso", &muonsRelIso);
t->SetBranchAddress("nMuons_CUT", &nMuons_CUT);
t->SetBranchAddress("nElectrons_CUT", &nElectrons_CUT);
t->SetBranchAddress("nJets_CUT", &nJets_CUT);
t->SetBranchAddress("ht", &ht);
t->SetBranchAddress("mht", &mht);
t->SetBranchAddress("met", &met);
t->SetBranchAddress("metphi", &metphi);
t->SetBranchAddress("dPhi0_CUT", &dPhi0_CUT);
t->SetBranchAddress("dPhi1_CUT", &dPhi1_CUT);
t->SetBranchAddress("dPhi2_CUT", &dPhi2_CUT);
 
int hadtauflag;
int nelectron,nmuon,njet;
double  dphi1 , dphi0 , dphi2;
TFile f1("selected.root","RECREATE");
TTree *nt = new TTree("model","model info");
TH1D *hjetspt=new TH1D("hjetspt", "jetspt histogram", 100,0,1000);
TH1D *hjetseta=new TH1D("hjetseta", "jetseta histogram", 10,-5.5,5.5);
TH1D *hjetsphi=new TH1D("hjetsphi", "jetsphi histogram", 10,-4,4);
TH1D *hgenmupt=new TH1D("hgenmupt", "genmupt histogram", 100, 0,100);
TH1D *hgenmueta=new TH1D("hgenmueta", "genmueta histogram", 10,-6,6);
TH1D *hgenmuphi=new TH1D("hgenmuphi", "genmuphi histogram", 10,-3.5,3.5);
TH1I *helec = new TH1I("helec","nelectron histogram",10,0,10);
TH1I *hjets =new TH1I("hjets","njets histogram",10,0,10);
TH1F *hht =new TH1F("hht","ht histogram",10,0,500);
TH1D *hmht =new TH1D("hmht","mht histogram",10,0,500);
TH1I *hmuon = new TH1I("hmuon","nmuon histogram",10,0,10);
TH1D *hmet =new TH1D("hmet","met histogram",10,0,500);
TH1D *hmetphi =new TH1D("hmetphi","nmetphi histogram",10,0,500);
TH1D *hdphi0 =new TH1D("hdphi0","dPhi0 histogram",10,0,4);
TH1D *hdphi1 =new TH1D("hdphi1","dPhi1 histogram",10,0,4);
TH1D *hdphi2 =new TH1D("hdphi0","dPhi2 histogram",10,0,4);
Int_t nentries =(Int_t)t->GetEntries();
cout<<nentries<<endl;

for(int j=0;j<nentries;j++){
t->GetEntry(j);
cout<<"jet size"<<jetsLVec->size()<<endl;
/*for (unsigned k=0;k<jetsLVec->size();k++){
hjetspt->Fill(jetsLVec[k].Pt());
hjetseta->Fill(jetsLVec[k].Eta());
hjetsphi->Fill(jetsLVec[k].Phi());
}*/
cout<<"genDecayStr : "<<genDecayStrVec->front().c_str()<<endl;
           for(int iv=0; iv<(int)genDecayLVec->size(); iv++){
              int pdgId = genDecayPdgIdVec->at(iv);
              printf("PDGId: %d,Index: %d/ Mom Index: %d  \n", pdgId, genDecayIdxVec->at(iv), genDecayMomIdxVec->at(iv));
            }

int NJET=0;
int genmu=0,accmu=0;
bool doneacc= false;
/*for (unsigned ij=0;ij<hjetspt->size();ij++){
if(hjetspt->at(ij)>30 && fabs(hjetseta->at(ij)<2.4))NJET++;
}*/
if(W_emuVec->size()!=0 && W_tau_emuVec->size()==0){

/*for(int l=0;l<genDecayLVec->size();l++){
int pdgId = genDecayPdgIdVec->at(l);
if(abs(pdgId)==13){     
hgenmupt->Fill(genDecayLVec[l].Pt());
hgenmueta->Fill(genDecayLVec[l].Eta());
hgenmuphi->Fill(genDecayLVec[l].Phi());
}
      if(abs(pdgId)==13){
        genmu++;
        double pt6= genDecayLVec->at(l).Pt();
        double eta6 = genDecayLVec->at(l).Eta();
        double phi6 = genDecayLVec->at(l).Phi();
        if(pt6>50 && fabs(eta6)<2.4)accmu++;
            doneacc=true;
        
       }
    }*/
}
 hjets->Fill(nJets_CUT);
 hht->Fill(ht);
 hmet->Fill(met);
 hmht->Fill(mht);
 hmetphi->Fill(metphi);
 hdphi0->Fill(dPhi0_CUT);
 hdphi1->Fill(dPhi1_CUT);
 hdphi2->Fill(dPhi2_CUT);

 nt->Fill();
}
f1.Write();

return 0;
}