#include "TStyle.h"
#include "TList.h"
#include "TH3.h"
#include "TH2.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 "TString.h"

#include <cstdlib>
#include<string> 
#include<iostream> 
#include<fstream>
#include<iomanip>
#include<cmath>
#include <vector>
#include <unistd.h>

void Th3DTry(){

  TChain *ch1 = new TChain("Rho");
  ch1->Add("run_56515_skimmed_tripinc_withCuts_0.root");
  Int_t entry = ch1->GetEntries();

  hprof2dpip = new TProfile2D("hprof2dpip","Profile of #pi^{+} Momenta versus Mass and BeamE", 100,0.45,0.55,100,0.0,2.8,1.1,5.7);
  hprof2dpim = new TProfile2D("hprof2dpim","Profile of #pi^{-} Momenta versus Mass and BeamE", 100,0.45,0.55,100,0.0,2.8,1.1,5.7);
  TH2D *h2d = new TH2D("h2d","h2d",100,0.45,0.55,100,0.0,2.8);

  Double_t IV_PipPim, DOCA_ang, E_g, Pip_Ptot, Pim_Ptot;
  ch1->SetBranchAddress("IV_PipPim",&IV_PipPim);  
  ch1->SetBranchAddress("DOCA_ang",&DOCA_ang);  
  ch1->SetBranchAddress("E_g",&E_g);
  ch1->SetBranchAddress("Pip_Ptot",&Pip_Ptot);
  ch1->SetBranchAddress("Pim_Ptot",&Pim_Ptot);

  for(int i =0;i<entry;i++){
    if(!(i%20000)) std::cout << "done " << i << " out of " << entry << " ==> " << double(i)*100.0/double(entry) << "%" << std::endl;  

    ch1->GetEntry(i);
    hprof2dpip->Fill(IV_PipPim,Pip_Ptot,E_g);
    hprof2dpim->Fill(IV_PipPim,Pim_Ptot,E_g);
    h2d->Fill(IV_PipPim,Pim_Ptot);


    
  }
  TCanvas *c1 = new TCanvas("c1","c1");
  c1->cd();
  gStyle->SetPalette(1);
  hprof2dpip->Draw("surf2");
  hprof2dpip->GetXaxis()->SetTitle("M(#pi^{+}#pi^{-}) [GeV]");
  hprof2dpip->GetXaxis()->SetTitleOffset(1.6);
  hprof2dpip->GetZaxis()->SetTitle("Beam Energy [GeV]");
  hprof2dpip->GetZaxis()->SetTitleOffset(1.6);
  hprof2dpip->GetYaxis()->SetTitle("#pi^{+} Momenta [GeV]");
  hprof2dpip->GetYaxis()->SetTitleOffset(1.3);

  TCanvas *c2 = new TCanvas("c2","c2");
  c2->cd();
  hprof2dpim->Draw("surf2");
  hprof2dpim->GetXaxis()->SetTitle("M(#pi^{+}#pi^{-}) [GeV]");
  hprof2dpim->GetXaxis()->SetTitleOffset(1.6);
  hprof2dpim->GetZaxis()->SetTitle("Beam Energy [GeV]");
  hprof2dpim->GetZaxis()->SetTitleOffset(1.6);
  hprof2dpim->GetYaxis()->SetTitle("#pi^{-} Momenta [GeV]");
  hprof2dpim->GetYaxis()->SetTitleOffset(1.3);

  TCanvas *c3 = new TCanvas("c3","c3");
  c3->cd();
  h2d->Draw("surf2");
  h2d->GetXaxis()->SetTitle("M(#pi^{+}#pi^{-}) [GeV]");
  h2d->GetXaxis()->SetTitleOffset(1.6);
  h2d->GetYaxis()->SetTitle("#pi^{+} Momenta [GeV]");
  h2d->GetYaxis()->SetTitleOffset(1.3);

  
}