#if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
#include <TSystem.h>
#include <TStyle.h>
#include <TFile.h>
#include <TTree.h>
#include <TCanvas.h>
#include <TH1D.h>
#include <TLorentzVector.h>
#include "TFile.h"
#include "TTree.h"
#include "TH2D.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TMath.h"
#include <iostream>
#endif

using namespace std;
#define PI 3.14159265359

namespace {
  int nfiles = 1;
  const char* inputs [] = {
       "trk_aug29_8M.root"  //input file
    
  };

  
}

void asymmetry() {
        gStyle->SetOptFit();
	TCanvas *c0 = new TCanvas("c0","c0"); c0->SetGrid();
	TCanvas *c1 = new TCanvas("c1","c1"); c1->SetGrid();
	
	TCanvas *c3 = new TCanvas("c3","c3"); c3->SetGrid();
        TCanvas *c4 = new TCanvas("c4","c4"); c4->SetGrid();
	TCanvas *c5 = new TCanvas("c5","c5"); c5->SetGrid();
	TCanvas *c6 = new TCanvas("c6","c6"); c6->SetGrid();
	TCanvas *c7 = new TCanvas("c7","c7"); c7->SetGrid();
        TCanvas *c8 = new TCanvas("c8","c8"); c8->SetGrid();
        TCanvas *c9 = new TCanvas("c9","c9"); c9->SetGrid();
       
	
	for (int i=0; i<nfiles; ++i) {
		 TH1D *hx = new TH1D("hx","hx",20,0,1); //Ndimu vs Xtarget
		 TH1D *hnum = new TH1D("hnum","hnum",10, 0,10);//Ndimu vs mass
                 TH1D *hphi_s = new TH1D("hphi_s","hphi_s",30, -1*PI,PI);//Ndimu vs phi_target_rest_frame
                 
                 TH1D *asym = new TH1D("Asym","Asym",5,0,0.5);//Sivers asymetry An

                //Ndimu vs phi_target_rest_frame for each bin of xtarget:
		 TH1D *x_target[10]; 
                 char name1[100];
                 for (int ii=0; ii<10; ii++) 
                  {
                    float xlow = 0.1*ii;
                    sprintf(name1,"xTarget_%f",xlow+0.05);
                    x_target[ii] = new TH1D(name1,name1, 20, -1*PI, PI);
                    x_target[ii]->Sumw2();
                  }

		//
		// Loop a TTree in a file
		
		TFile *f = TFile::Open(inputs[i],"read");
		TTree *T = (TTree*) f->Get("Truth");

#define MAX_PARTICLES 100
		int gnhodo[MAX_PARTICLES];
		float gpx[MAX_PARTICLES], gpy[MAX_PARTICLES], gpz[MAX_PARTICLES];
		int n_particles, krecstat;

		T->SetBranchAddress("gnhodo", &gnhodo);
		T->SetBranchAddress("gpx", &gpx);
		T->SetBranchAddress("gpy", &gpy);
		T->SetBranchAddress("gpz", &gpz);
		
		T->SetBranchAddress("krecstat", &krecstat);

		float mu_mass = .105658;
                float pro_mass = 0.93827;
                float beamE = 120.0;
                float gE[2];
                float phi, phi_s;
		
		for(int ientry=0; ientry<T->GetEntries(); ++ientry) {
			T->GetEntry(ientry);
                        
			//Initial 4-vector definition
			TLorentzVector mu0, mu1, dimu, hadron_b, hadron_t, initial,spin;
                        gE[0] = TMath::Sqrt(gpx[0]*gpx[0]+gpy[0]*gpy[0]+gpz[0]*gpz[0]+mu_mass*mu_mass);
                        gE[1] = TMath::Sqrt(gpx[1]*gpx[1]+gpy[1]*gpy[1]+gpz[1]*gpz[1]+mu_mass*mu_mass);
			mu0.SetPxPyPzE(gpx[0], gpy[0], gpz[0], gE[0]);
			mu1.SetPxPyPzE(gpx[1], gpy[1], gpz[1], gE[1]);
                        hadron_b.SetPxPyPzE(0.0,0.0,TMath::Sqrt(beamE*beamE-pro_mass*pro_mass),beamE);
                        hadron_t.SetPxPyPzE(0.0,0.0,0.0,pro_mass);
         		if (ientry % 2 ==0)
                       {
                        spin.SetPxPyPzE(0.0,1.0,0.0,0.0);
                       }
             		else
                        {
                        spin.SetPxPyPzE(0.0,-1.0,0.0,0.0);
                        }
                        initial = hadron_b + hadron_t;
                        dimu = mu0+mu1;
			float dimu_gmass = dimu.Mag();
                      
                        
                      //boost to centre of mass

                        TLorentzRotation cm_r;
                        TVector3 cm_v = initial.BoostVector();
                        cm_r.Boost(-cm_v);

                        TLorentzVector hadron_b_cm = cm_r * hadron_b;
                        TLorentzVector hadron_t_cm = cm_r * hadron_t;
                        TLorentzVector mu0_cm = cm_r * mu0;
                        TLorentzVector mu1_cm = cm_r * mu1;
                        TLorentzVector dimu_cm = cm_r * dimu;
                        TLorentzVector spin_cm = cm_r * spin;

                        float s = 2.0*beamE*pro_mass + 2.0*pro_mass*pro_mass;
                        float sqrt_s = TMath::Sqrt(s);
                        float dimu_E = dimu_cm.E();
                        float dimu_pl = dimu_cm.Pz();
                        float dimu_pl_max = sqrt_s/2.0*(1-(dimu_gmass*dimu_gmass)/s);
                        float xf = dimu_pl/dimu_pl_max;
                        float y = 0.5*TMath::Log((dimu_E+dimu_pl)/(dimu_E-dimu_pl));
                        float tau = (dimu_gmass*dimu_gmass)/s;
                        float xb = 0.5*(TMath::Sqrt(xf*xf+4*tau)+xf);
                        float xt = 0.5*(TMath::Sqrt(xf*xf+4*tau)-xf);
                        
                        //boost to dilepton restframe

                        TLorentzRotation cm_dimu;
                        TVector3 cm_v2 = dimu_cm.BoostVector();
                        cm_dimu.Boost(-cm_v2);

                        TLorentzVector hadron_b_cm_dimu = cm_dimu * hadron_b_cm;
                        TLorentzVector hadron_t_cm_dimu = cm_dimu * hadron_t_cm;
                        TLorentzVector mu0_cm_dimu = cm_dimu * mu0_cm;
                        TLorentzVector mu1_cm_dimu = cm_dimu * mu1_cm;
                        TLorentzVector dimu_cm_dimu = cm_dimu * dimu_cm;
                        TLorentzVector spin_cm_dimu = cm_dimu * spin_cm;


                        TVector3 zaxis = (hadron_b_cm_dimu.Vect()).Unit()-(hadron_t_cm_dimu.Vect()).Unit();
                        TVector3 zunit = zaxis.Unit();
                        TVector3 hadron_plane_norm = zunit.Cross((hadron_b_cm_dimu.Vect()).Unit());
                        TVector3 lepton_plane_norm = zunit.Cross((mu0_cm_dimu.Vect()).Unit());
                        TVector3 yunit = lepton_plane_norm.Unit();
                        TVector3 xunit = (yunit.Cross(zunit)).Unit();
                 
                                           
                        

                        
                     //below are phi_s_tf from target rest frame. The frame for the Asymmetry calculation
                         TVector3 spin_vect_tf = spin.Vect();
                         TVector3 z_prime_unit = (hadron_b.Vect()).Unit();
                         TVector3 y_prime_unit = (z_prime_unit.Cross(dimu.Vect())).Unit(); 
                         TVector3 x_prime_unit = (y_prime_unit.Cross(z_prime_unit)).Unit(); 
                         float phi_s_tf;                     
                         float cosphi_s_tf = (x_prime_unit.Dot(spin_vect_tf))/((x_prime_unit.Mag())*(spin_vect_tf.Mag()));                 
                         if ( acos( spin_vect_tf.Dot(y_prime_unit) / ( (spin_vect_tf).Mag() * y_prime_unit.Mag() ) ) > (PI/2) )       
                           {
                            phi_s_tf = -acos(cosphi_s_tf);
                          }
                        else
                          {
                            phi_s_tf = acos(cosphi_s_tf);
                          } 

                          int bin_xtarget = xt/0.1;
                 cerr<<ientry<<"  "<<T->GetEntries()<<" "<<dimu_gmass<<" "<<xt<<" "<<phi_s_tf<<" "<<endl;       

//good_event status -> currently only for successfully reconstructed event. We can also add trigger condition later for the good event status:
	       bool good_event=false;
               if (krecstat==0)
                {
                  good_event=true; 
                }  
                           
                        if(good_event) //only for successfully reconstructed event
                         {
			   hnum->Fill(dimu_gmass);
                           hx->Fill(xt);
                           hphi_s->Fill(phi_s_tf);
                           
                          if (bin_xtarget >=0 && bin_xtarget < 10)
                            {
                             x_target[bin_xtarget]->Fill(phi_s_tf);
                               
                            }
                         }
		}

		// End of loop
                cerr<<"end of loop"<<endl;
                cerr<<endl;
		cerr<<"get An"<<endl;
               //Asymmetry calculation:
                for (int j=0; j<5; j++)
                 {
                   float sum=0;
                   float sum_N=0;
                   float An=0;
                   float num_err=0;
                   float den_err=0;
                   float tot_err=0;
    	 	   float new_err=0;
                   
                   for (int k=1; k<=20; k++)
                     {
                      float phi_angle = -1.0*PI+(k-1)*2.0*PI/20.0+0.5*2.0*PI/20.0;
                      float a = x_target[j]->GetBinContent(k);
                      float b = x_target[j]->GetBinError(k);
                      sum = sum + a*sin(phi_angle);
                      sum_N = sum_N+a;
                      num_err = num_err + (b*sin(phi_angle))*(b*sin(phi_angle));
                      den_err = den_err + b*b;
		      //New error calculation from Kenichi. Kenichi pointed out that they are correlated:
		      new_err = new_err + pow(b * (sin(phi_angle) * sum_N - sum) / pow(sum_N, 2), 2);
                      
                     }
                    num_err = TMath::Sqrt(num_err);
                    den_err = TMath::Sqrt(den_err);
                    An = 2*sum/sum_N;
                    tot_err = 2*An*TMath::Sqrt((num_err/sum)*(num_err/sum)+(den_err/sum_N)*(den_err/sum_N));
                    
		   if(fabs(An)<=1)
                    {
                     asym->SetBinContent(j+1,An);
                     //asym->SetBinError(j+1,tot_err);
                     //applying new error calculation from kenichi:
                     new_err = 2*TMath::Sqrt(new_err);
                     asym->SetBinError(j+1,new_err);
                    }
		    cerr<<"Bin Xtarget = "<<j+1<<"  An = "<<An<<endl;
                 }
                     

		//some cosmetics
		hnum->SetMarkerStyle(20);
		hx->SetMarkerStyle(20);
		hphi_s->SetMarkerStyle(20);
                asym->SetMarkerStyle(20);
                for (int j=0; j<10; j++)
                 {
                   x_target[j]->SetMarkerStyle(20);
                 }
		

		int color = kBlack;

		if(i==0) {

			color = kBlack;

			c0->cd();
			//c0->SetLogy();
			hnum->SetTitle("nDimu-trig. vs. M_{#mu#mu}; M_{#mu#mu} [GeV/c^{2}]; nDimu-trig.");
			hnum->SetMarkerColor(color);
			hnum->SetLineColor(color);
			hnum->SetMinimum(1);
			hnum->Draw("e");
			hnum->SetStats(0);

			c1->cd();
			//c1->SetLogy();
			hx->SetTitle("nDimu-trig. vs. x_target; x_target; nDimu-trig.");
			hx->SetMarkerColor(color);
			hx->SetLineColor(color);
			hx->SetMinimum(1);
			hx->Draw("e");
			hx->SetStats(0);

			

                        c3->cd();
			hphi_s->SetTitle("nDimu-trig. vs. #phi_{s}; #phi_{s}; nDimu-trig.");
			hphi_s->SetMarkerColor(color);
			hphi_s->SetLineColor(color);
			hphi_s->SetMinimum(0.);
			//hphi->SetMaximum(1.1);
			hphi_s->Draw("e");
			hphi_s->SetStats(0);

                        c4->cd();
			x_target[0]->SetTitle("nDimu-trig. vs. #phi_{s} (xt = 0.05) ; #phi_{s}; nDimu-trig.");
			x_target[0]->SetMarkerColor(color);
			x_target[0]->SetLineColor(color);
			x_target[0]->SetMinimum(0.);
			x_target[0]->Draw("e");
                        x_target[0]->SetStats(0);

                        c5->cd();
			x_target[1]->SetTitle("nDimu-trig. vs. #phi_{s} (xt = 0.15) ; #phi_{s}; nDimu-trig.");
			x_target[1]->SetMarkerColor(color);
			x_target[1]->SetLineColor(color);
			x_target[1]->SetMinimum(0.);
			x_target[1]->Draw("e");
			x_target[1]->SetStats(0);

                        c6->cd();
			x_target[2]->SetTitle("nDimu-trig. vs. #phi_{s} (xt = 0.25) ; #phi_{s}; nDimu-trig.");
			x_target[2]->SetMarkerColor(color);
			x_target[2]->SetLineColor(color);
			x_target[2]->SetMinimum(0.);
			x_target[2]->Draw("e");
			x_target[2]->SetStats(0);

                        c7->cd();
			x_target[3]->SetTitle("nDimu-trig. vs. #phi_{s} (xt = 0.35) ; #phi_{s}; nDimu-trig.");
			x_target[3]->SetMarkerColor(color);
			x_target[3]->SetLineColor(color);
			x_target[3]->SetMinimum(0.);
			x_target[3]->Draw("e");
			x_target[3]->SetStats(0);

                        c8->cd();
			x_target[4]->SetTitle("nDimu-trig. vs. #phi_{s} (xt = 0.45) ; #phi_{s}; nDimu-trig.");
			x_target[4]->SetMarkerColor(color);
			x_target[4]->SetLineColor(color);
			x_target[4]->SetMinimum(0.);
			x_target[4]->Draw("e");
			x_target[4]->SetStats(0);

                        c9->cd();
			asym->SetTitle("Drell-Yan Target Single-Spin Asymmetry; X_{target}; A_{N}");
			asym->SetMarkerColor(color);
			asym->SetLineColor(color);
			asym->SetMinimum(-0.5);
			asym->SetMaximum(0.5);
			asym->Draw("e");
			asym->SetStats(0);

                       

			
		}
		else {
			color = i+1;

			c0->cd();
			hnum->SetMarkerColor(color);
			hnum->SetLineColor(color);
			hnum->Draw("esame");

			c1->cd();
			hx->SetMarkerColor(color);
			hx->SetLineColor(color);
			hx->Draw("esame");

			

                        c3->cd();
			hphi_s->SetMarkerColor(color);
			hphi_s->SetLineColor(color);
			hphi_s->Draw("esame");

			
		}
	}


}