#define WeightAnalysis_Res__ntuple_cxx
#include "WeightAnalysis_Res__ntuple.h"
#include <TH2.h>
#include <TStyle.h>
#include <TCanvas.h>

void WeightAnalysis_Res__ntuple::Loop()
{
//   In a ROOT session, you can do:
//      root> .L WeightAnalysis_Res__ntuple.C
//      root> WeightAnalysis_Res__ntuple t
//      root> t.GetEntry(12); // Fill t data members with entry number 12
//      root> t.Show();       // Show values of entry 12
//      root> t.Show(16);     // Read and show values of entry 16
//      root> t.Loop();       // Loop on all entries
//

//     This is the loop skeleton where:
//    jentry is the global entry number in the chain
//    ientry is the entry number in the current Tree
//  Note that the argument to GetEntry must be:
//    jentry for TChain::GetEntry
//    ientry for TTree::GetEntry and TBranch::GetEntry
//
//       To read only selected branches, Insert statements like:
// METHOD1:
//    fChain->SetBranchStatus("*",0);  // disable all branches
//    fChain->SetBranchStatus("branchname",1);  // activate branchname
// METHOD2: replace line
//    fChain->GetEntry(jentry);       //read all branches
//by  b_branchname->GetEntry(ientry); //read only this branch
   if (fChain == 0) return;

   Long64_t nentries = fChain->GetEntriesFast();

   Long64_t nbytes = 0, nb = 0;
   for (Long64_t jentry=0; jentry<nentries;jentry++) {
      Long64_t ientry = LoadTree(jentry);
      if (ientry < 0) break;
      nb = fChain->GetEntry(jentry);   nbytes += nb;
      // if (Cut(ientry) < 0) continue;
      //std::cout << " Running " << std::endl;
      
      int n_muons = muonsE->size();
      int n_jets = jetsE->size();
      double leadingJetPtCut = 30;
      double metCut = 55;
      double met_HPTOCut = 55;
      double x1Cut_1 = 0.1;
      double x1Cut_2 = 1.0;
      double x2Cut_1 = 0.1;
      double x2Cut_2 = 1.0;
      bool useHPTOmet = true;
      double met,pmissx,pmissy;

      vector<TLorentzVector> finalMuons;
      for (int i=0;i<n_muons;i++){
         TLorentzVector xMuon;
         xMuon.SetPtEtaPhiE(muonsPt->at(i),muonsEta->at(i),muonsPhi->at(i),muonsE->at(i));
         finalMuons.push_back(xMuon);
      }

      vector<TLorentzVector> finalJets;
      for (int i=0;i<n_jets;i++){
         TLorentzVector xJet;
         xJet.SetPtEtaPhiE(jetsPt->at(i),jetsEta->at(i),jetsPhi->at(i),jetsE->at(i));
         finalJets.push_back(xJet);
      }

      vector<TLorentzVector> finalTruthJets;
      for (int i=0;i<n_jets;i++){
         TLorentzVector xJetTruth;
         xJetTruth.SetPtEtaPhiE(jetsTruthPt->at(i),jetsTruthEta->at(i),jetsTruthPhi->at(i),jetsTruthE->at(i));
         finalTruthJets.push_back(xJetTruth);
      }

      double ptMuon1,ptMuon2,met_HPTOx,met_HPTOy,met_HPTO;
      double x1,x2,collmass;
      if (n_muons > 1) {
        ptMuon1 = finalMuons.at(0).Pt();
        ptMuon2 = finalMuons.at(1).Pt();
        // MET HPTO
        met_HPTOx= - finalMuons.at(0).Px() - finalMuons.at(1).Px();
        met_HPTOy= - finalMuons.at(0).Py() - finalMuons.at(1).Py();
        for (int i = 0 ; i < finalJets.size() ; i++){
            met_HPTOx += - finalJets.at(i).Px();
            met_HPTOy += - finalJets.at(i).Py();
        }
        met_HPTO = sqrt(met_HPTOx*met_HPTOx + met_HPTOy*met_HPTOy);
      
        met=met_HPTO;
        pmissx = met_HPTOx;
        pmissy = met_HPTOy;
        double collmas;
        // COLLINEAR APPROXIMATION
        if (useHPTOmet){
            std::cout << " Point 1 " << std::endl;
            std::cout << " finalMuons.at(0).Pt() " << finalMuons.at(0).Pt() << std::endl;
            std::cout << " finalMuons.at(1).Pt() " << finalMuons.at(1).Pt() << std::endl;
            std::cout << " met_HPTOx " << met_HPTOx << std::endl;
            std::cout << " met_HPTOy " << met_HPTOy << std::endl;

            resultMass = MassCollinearCore(finalMuons.at(0), finalMuons.at(1), //particles
                                    met_HPTOx, met_HPTOy,
                                    collmass, x1, x2);
            std::cout << " Point 2 " << std::endl;
        } else {
            //resultMass = MassCollinearCore(finalMuons.at(0), finalMuons.at(1), //particles
            //                        pmissx, pmissy,
            //                        collmass, x1, x2);
        }

      }  


      double weightHistograms = eventWeight;
      WeightAnalysis_Res__eventWeight_0->Fill(eventWeight,weightHistograms);
      WeightAnalysis_Res__n_muons_0->Fill(n_muons,weightHistograms);

      if (!(muonsE->size()==2 && muonsCharge->at(0)*muonsCharge->at(1)< 0.)) continue;

      WeightAnalysis_Res__n_muons_1->Fill(n_muons,weightHistograms);

      if (!((finalMuons.at(0)+finalMuons.at(1)).M()>30. && (finalMuons.at(0)+finalMuons.at(1)).M()< 75.)) continue;

      WeightAnalysis_Res__n_muons_2->Fill(n_muons,weightHistograms);

      if (!(n_jets > 0)) continue;

      WeightAnalysis_Res__n_muons_3->Fill(n_muons,weightHistograms);

      if (!(jetsPt->at(0) > leadingJetPtCut)) continue;

      WeightAnalysis_Res__n_muons_4->Fill(n_muons,weightHistograms);

      if (!(met > metCut)) continue;

      WeightAnalysis_Res__n_muons_5->Fill(n_muons,weightHistograms);

      if (!(met_HPTO > met_HPTOCut)) continue;

      WeightAnalysis_Res__n_muons_6->Fill(n_muons,weightHistograms);

      if (!(( x1 > x1Cut_1 &&  x1 < x1Cut_2) && (x2 > x2Cut_1 && x2 < x2Cut_2 ))) continue;

      WeightAnalysis_Res__n_muons_7->Fill(n_muons,weightHistograms);

   }
   WeightAnalysis_Res__n_muons_0->Draw();
}