#include <TFile.h>
#include <TTree.h>
#include <TROOT.h>
#include <TH1.h>
//#include "UserCode/TopAnalysis/interface/MiniEvent.h"
//#include "UserCode/TopAnalysis/interface/ReadTree.h"
#include <TH2.h>
#include <TSystem.h>
#include <vector>
#include <TLorentzVector.h>
#include <math.h>
#include "TMath.h"
#include <iostream>
#include <fstream>
#include <algorithm>
#include <map>
#include <TFileMerger.h>

void lorentz32()
//void stupid_TLV(TString filename,  TString outname)

{
   TFile *f = new TFile("May26_2.root");
   TTree *t = (TTree*)f->Get("demo/AnaTree");
   TH1F *h1 = new TH1F("h1","Z1_mass",100,0,300);
	TH1F *h2 = new TH1F("h2","Z2_mass",100,0,300);
	TH1F *h3 = new TH1F("h3","Z3_mass",100,0,300);

using namespace std;

float SelMu_px, SelMu_py, SelMu_pz, SelMu_en, SelMu_pt, SelMu_eta, SelMu_phi;
int SelMu, SelMu_charge;

          t->SetBranchAddress("SelMu_px",&SelMu_px);
          t->SetBranchAddress("SelMu_py",&SelMu_py);
          t->SetBranchAddress("SelMu_pz",&SelMu_pz);
          t->SetBranchAddress("SelMu_en",&SelMu_en);
          t->SetBranchAddress("SelMu",&SelMu);
          t->SetBranchAddress("SelMu_charge",&SelMu_charge);

 TLorentzVector tmp(0,0,0,0), tmp1(0,0,0,0), tmp2(0,0,0,0), tmp3(0,0,0,0);
 for (int i=0;i<10000;i++)

           {
 TLorentzVector tmp, tmp1, tmp2, tmp3, tmp4;
         tmp.SetPxPyPzE(SelMu_px,SelMu_py,SelMu_pz,SelMu_en);
         t->GetEntry(i);

	 cout<< " i : "<<i<<" nmu : "<<SelMu<<endl;
    
 for(int m=0; m<SelMu;m++)
        {

	  if(SelMu<2)continue;
  
int m1=0, m2=1, m3=2;

// tmp1.SetPxPyPzE(SelMu_px[m0], SelMu_py[m0], SelMu_pz[m0], SelMu_en[m0]);
if(m==0) tmp1.SetPxPyPzE(SelMu_px[m], SelMu_py[m], SelMu_pz[m], SelMu_en[m]);
if(m==1) tmp2.SetPxPyPzE(SelMu_px[m], SelMu_py[m], SelMu_pz[m], SelMu_en[m]);
if(m==2) tmp3.SetPxPyPzE(SelMu_px[m], SelMu_py[m], SelMu_pz[m], SelMu_en[m]);


//if (SelMu_charge[m2]*SelMu_charge[m3]>0.) continue;
//if (SelMu_charge[m3]*SelMu_charge[m1]<0.); continue;

if(m<1 )continue;
if (SelMu_charge[1]*SelMu_charge[2]>0.); continue;
 
//      double z_mass = (tmp1 + tmp2).Pt();  
//      double z_mass = (tmp2 + tmp3).Pt();

	double z1_mass = (tmp1 + tmp2).M();
	double z2_mass = (tmp2 + tmp3).M();
	double z3_mass = (tmp1 + tmp3).M();
	cout<< z1_mass<<endl;
   cout<< z2_mass<<endl;
   cout<< z3_mass<<endl;

//      double z_mass = z1_mass.M(); 
//      float z_mass = (z1_mass).M();
std::cout<< " z1_mass : " <<z1_mass<<" px : "<< SelMu_px[m1]<<  " py : "<< SelMu_py[m1] <<" "<< SelMu_pz[m1]<<"   "<< SelMu_en[m]<<std::endl;
   std::cout<< " z2_mass : " <<z1_mass<<" px : "<< SelMu_px[m2]<<  " py : "<< SelMu_py[m2] <<" "<< SelMu_pz[m2]<<std::endl;
   std::cout<< " z3_mass : " <<z1_mass<<" px : "<< SelMu_px[m3]<<  " py : "<< SelMu_py[m3] <<" "<< SelMu_pz[m3]<<std::endl;
std::cout<< " z2_mass : " <<z2_mass<<std::endl;
std::cout<< " z3_mass : " <<z3_mass<<std::endl;

         t->Fill();
         h1->Fill(z1_mass);
         h2->Fill(z2_mass);
         h3->Fill(SelMu_en[m]);
        
}
        }       
         TCanvas *c1 = new TCanvas("c1","c1",800,600) ;
        c1->Divide(2,2); 
  	c1->cd(1);
         h1->Draw();
         h1->GetXaxis()->SetTitle("Z_12_mass");
	 c1->cd(2);
         h2->Draw();
         h2->GetXaxis()->SetTitle("Z_23_mass");
  	c1->cd(3);
         h3->Draw();
         h3->GetXaxis()->SetTitle("Z_31_mass");

//         f->Write();
}