#include "TROOT.h"
#include "TRint.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TF1.h"
#include "TCanvas.h"
#include "TChain.h"
#include "TVector3.h"
#include "TMath.h"
#include "TFile.h"
#include "TSystem.h"
#include "TGStatusBar.h"
#include "TSystem.h"
#include "TXMLEngine.h"
#include "TTree.h"
#include "TLorentzVector.h"
#include "TNtuple.h" 
#include <cstdlib>
#include<string> 
#include<iostream> 
#include<fstream>
#include<iomanip>
#include<cmath>
#include <vector>
#include "TStyle.h"
#include "TGraphErrors.h"
#include "TGraph.h"
#include "TExec.h"
void Up_Stream(){
  Double_t M_P = 0.938272;   //Proton

  //############################################################################################################################################  
  //Setting options for histogram display
  gStyle->SetOptStat("ne");
  gStyle->SetOptFit(111);
  //############################################################################################################################################  
  //Calling data 
  TChain *chain = new TChain("Rho");  
  //chain->Add("/Volumes/Mac_Storage/YIELD/SINGLE_TRACK/YIELD_DATA/PRODUCTION/RUN_56515/Anal_Upstream0.root");
 chain->Add("small.root");

  //############################################################################################################################################  
  //HISTOGRAMS  
  Int_t nEvent = chain->GetEntries();
  //############################################################################################################################################   
  //Setting up Histograms to be used for profiling the beam correction
  //-------------Bin numbers and range are defined here
  Double_t pt_xlo = 0.0;
  Double_t pt_xup = 72;
  const Int_t pt_xnbin = 12;
  //Y Bins
  Double_t pt_ylo = -180;
  Double_t pt_yup = 180.;
  const Int_t pt_ynbin = 24;
  //Z Bins
  Double_t pt_zlo = 0.8;
  Double_t pt_zup = 1.2;
  const Int_t pt_znbin = 100;
  TH3D *PhiThetaMassPip = new TH3D("PhiThetaMassPip","PhiThetaMassPip",pt_xnbin,pt_xlo,pt_xup,pt_ynbin,pt_ylo,pt_yup,pt_znbin,pt_zlo,pt_zup);
  TH3D *PhiThetaMassPim = new TH3D("PhiThetaMassPim","PhiThetaMassPim",pt_xnbin,pt_xlo,pt_xup,pt_ynbin,pt_ylo,pt_yup,pt_znbin,pt_zlo,pt_zup);
  TH2D *PhiThetaPip = new TH2D("PhiThetaPip","PhiThetaPip",pt_xnbin,pt_xlo,pt_xup,pt_ynbin,pt_ylo,pt_yup);
  TH2D *PhiThetaPim = new TH2D("PhiThetaPim","PhiThetaPim",pt_xnbin,pt_xlo,pt_xup,pt_ynbin,pt_ylo,pt_yup);
  TH2D *PhiThetaPipDiff = new TH2D("PhiThetaPipDiff","PhiThetaPipDiff",pt_xnbin,pt_xlo,pt_xup,pt_ynbin,pt_ylo,pt_yup);
  TH2D *PhiThetaPimDiff = new TH2D("PhiThetaPimDiff","PhiThetaPimDiff",pt_xnbin,pt_xlo,pt_xup,pt_ynbin,pt_ylo,pt_yup);
  TH2D *PhiThetaPipSpect = new TH2D("PhiThetaPipSpect","PhiThetaPipSpect",pt_xnbin,pt_xlo,pt_xup,pt_ynbin,pt_ylo,pt_yup);
  TH2D *PhiThetaPimSpect = new TH2D("PhiThetaPimSpect","PhiThetaPimSpect",pt_xnbin,pt_xlo,pt_xup,pt_ynbin,pt_ylo,pt_yup);
  
  //############################################################################################################################################   
//READ-----------I TOOK OUT DATA TO USE EXCLUSIVIVTY BECAUSE I USED COPYTREE TO SKIM AN EXCLUSIVE DATA SET   

  Double_t Pip_Phi, Pip_Theta, Pim_Theta, Pim_Phi, mm_PipPim;
  chain->SetBranchAddress("Pip_Phi",&Pip_Phi);
  chain->SetBranchAddress("Pip_Theta",&Pip_Theta);
  chain->SetBranchAddress("Pim_Theta",&Pim_Theta);
  chain->SetBranchAddress("Pim_Phi",&Pim_Phi);
  chain->SetBranchAddress("mm_PipPim",&mm_PipPim);


  //############################################################################################################################################   
  //-------------Analyzing Data
  for( Int_t i = 0; i < nEvent ; i++){//nEvent
    chain->GetEntry(i); 
    if(!(i%100000)) std::cout << "done " << i << " out of " << nEvent << " ==> " << double(i)*100.0/double(nEvent) << "%" << std::endl;  
    //############################################################################################################################################   
    //-------------Applying cuts and filling histograms
      PhiThetaMassPip->Fill(Pip_Theta,Pip_Phi,mm_PipPim); 
      PhiThetaMassPim->Fill(Pim_Theta,Pim_Phi,mm_PipPim); 
      PhiThetaPipSpect->Fill(Pip_Theta,Pip_Phi);
      PhiThetaPimSpect->Fill(Pim_Theta,Pim_Phi);

  }
  //############################################################################################################################################   

  TH1D *Proj_[pt_xnbin][pt_ynbin]; //<========
  TH1D *Proj2_[pt_xnbin][pt_ynbin]; //<========

  //TCanvas *Can1_[pt_xnbin][pt_ynbin];
  //TCanvas *Can2_[pt_xnbin][pt_ynbin];

  for (Int_t cx = 0; cx < pt_xnbin; cx++) {
      for (Int_t cy = 0; cy < pt_ynbin ; cy++) {//pt_ynbin
        
        char c[10]; char c2[10];
        sprintf(c, "%d",cx);
        sprintf(c2, "%d",cy);
        TString nam0 = "Proj_";
        TString nam01 = "Proj2_";

        TString nam1 = "Can1_";
        TString nam2 = "Can2_";
        TString nam3 = nam0+c+c2;
        TString nam4 = nam1+c+c2;
        TString nam5 = nam01+c+c2;
        TString nam6 = nam2+c+c2;
        
        //Can1_[cx][cy]= new TCanvas(nam4,nam4);
        Proj_[cx][cy]=PhiThetaMassPip->ProjectionZ(nam3,cx,cx+1,cy,cy+1);
 
        //Can2_[cx][cy]= new TCanvas(nam6,nam6);
        Proj2_[cx][cy]=PhiThetaMassPim->ProjectionZ(nam5,cx,cx+1,cy,cy+1);
        
        Double_t initialPar = Proj_[cx][cy]->GetMean(); 
        Double_t initial2Par = Proj2_[cx][cy]->GetMean(); 

        
        TF1 fitter("fitter","gaus + [3] + [4]*x + [5]*x*x + [6]*x*x*x + [7]*x*x*x*x",0.9,1.0);
        fitter.SetParName(0, "Norm"); fitter.SetParName(1, "Mean"); fitter.SetParName(2, "#sigma"); fitter.SetParName(3, "const term");
        fitter.SetParName(4, "linear term"); fitter.SetParName(5, "quad term"); fitter.SetParName(6, "cube term");
        fitter.SetParameters(24000, initialPar, 0.01, -1.5, 0, 100, 100,100);
        fitter.SetParLimits(0,10,80000); fitter.SetParLimits(1,initialPar-0.05,initialPar+0.05);

        TF1 fitter2("fitter2","gaus + [3] + [4]*x + [5]*x*x + [6]*x*x*x + [7]*x*x*x*x",0.9,1.0);
        fitter2.SetParName(0, "Norm"); fitter2.SetParName(1, "Mean"); fitter2.SetParName(2, "#sigma"); fitter2.SetParName(3, "const term");
        fitter2.SetParName(4, "linear term"); fitter2.SetParName(5, "quad term"); fitter2.SetParName(6, "cube term");
        fitter2.SetParameters(24000, initial2Par, 0.01, -1.5, 0, 100, 100,100);
        fitter2.SetParLimits(0,10,80000); fitter2.SetParLimits(1,initial2Par-0.05,initial2Par+0.05);
        
        Proj_[cx][cy]->Fit("fitter","R,Q,N,0");  //USE Q for quiet ,N,0 
        Double_t temp1 = fitter.GetParameter(1);
        Double_t Para1;
        if(temp1 <0.7){Para1 = 0.0;}
        else{Para1 = M_P - temp1;}
        PhiThetaPip->SetBinContent(cx+1,cy+1,fitter.GetParameter(1));
        PhiThetaPip->SetBinError(cx+1,cy+1,fitter.GetParError(1));
        PhiThetaPipDiff->SetBinContent(cx+1,cy+1,Para1);
        PhiThetaPipDiff->SetBinError(cx+1,cy+1,fitter.GetParError(1));

        Proj2_[cx][cy]->Fit("fitter2","R,Q,N,0");  //USE Q for quiet ,N,0  
        Double_t temp2 = fitter2.GetParameter(1);
        Double_t Para2;
        if(temp2 <0.7){Para2 = 0.0;}
        else{Para2 = M_P - temp2;}
        PhiThetaPim->SetBinContent(cx+1,cy+1,fitter2.GetParameter(1));
        PhiThetaPim->SetBinError(cx+1,cy+1,fitter2.GetParError(1));
        PhiThetaPimDiff->SetBinContent(cx+1,cy+1,Para2);
        PhiThetaPimDiff->SetBinError(cx+1,cy+1,fitter2.GetParError(1));

      }
   }
  TCanvas *C1 = new TCanvas("C1","C1",800,650);
  TCanvas *C2 = new TCanvas("C2","C2",800,650);
  TExec *ex2 = new TExec("ex2","gStyle->SetPaintTextFormat(\".3f\");");
  TExec *ex1 = new TExec("ex1","gStyle->SetPaintTextFormat(\".0f\");");

  
//  PhiThetaPip->SetTitle("Mass of Proton Binned in #pi^{+} Angles");
//  PhiThetaPip->GetXaxis()->SetTitle("#pi^{+}#theta[#circ]");
//  PhiThetaPip->GetYaxis()->SetTitle("#pi^{+}#phi[#circ]");
  
  C1->Divide(2,1);C2->Divide(2,1);
  C1->cd(1);
  gPad->DrawFrame(pt_xlo,pt_ylo,pt_xup,pt_yup); ex2->Draw();
  PhiThetaPip->Draw("TEXT SAME");

  
  C1->cd(2);
  gPad->DrawFrame(pt_xlo,pt_ylo,pt_xup,pt_yup); ex1->Draw();
  gPad->SetLogz();
  PhiThetaPipSpect->Draw("TEXT SAME COLZ");
	//C1->Modified();
	C1->Update();
  
  C2->cd(1);
  gPad->DrawFrame(pt_xlo,pt_ylo,pt_xup,pt_yup); ex2->Draw();
  PhiThetaPim->Draw("TEXT SAME");

  C2->cd(2);
  gPad->DrawFrame(pt_xlo,pt_ylo,pt_xup,pt_yup); ex1->Draw();
  gPad->SetLogz();
  PhiThetaPimSpect->Draw("TEXT SAME COLZ");
  
  TCanvas *C3 = new TCanvas("C3","C3",800,650);
  TCanvas *C4 = new TCanvas("C4","C4",800,650);
  C3->cd();
  gPad->DrawFrame(pt_xlo,pt_ylo,pt_xup,pt_yup); ex2->Draw();
  PhiThetaPipDiff->Draw("TEXT SAME");
  C4->cd();
  gPad->DrawFrame(pt_xlo,pt_ylo,pt_xup,pt_yup); ex2->Draw();
  
  //##########Setting title#######
  PhiThetaPimDiff->SetTitle("Difference of Proton mass as Binned in #pi^{-} Angles");
  PhiThetaPimDiff->GetXaxis()->SetTitle("#pi^{-}#theta[#circ]");
  PhiThetaPimDiff->GetYaxis()->SetTitle("#pi^{-}#phi[#circ]");
  PhiThetaPimDiff->Draw("TEXT SAME");

  
}