//draw results of fast simulation analysis
//in nice and complete form
//16/04/2014 A.Karavdina

void DrawAnaResults(TString path = "/panda/scrut14/macro/scrut/fastsim_hc/", TString file_sig_name = "/hc_rho0pi0_fast.root_taskana.root", 
		    TString file_bkg_name ="/dpm_5_6119_fast.root_taskana.root", 
		    TString file_fs_name = "/hc_rho0pi0_5_6119_fast.root"){
 
  //  TString path = "/panda/scrut14/macro/scrut/fastsim_hc/";
  TString file_fs = path+file_fs_name;
  TString file_sig = path+file_sig_name;
  TString file_bkg = path+file_bkg_name;
  TFile *finFS = new TFile(file_fs.Data(),"READ");
  TTree *tree_final_state = (TTree*)finFS->Get("cbmsim");

  TH1 *hgammath = new TH1D("hgammath","#gamma ;#theta_{LAB}, ^#circ",1e2,0,180);
  tree_final_state->Project("hgammath","PidNeutralCand.GetMomentum().Theta()*TMath::RadToDeg()");
  TH1 *hpichth = new TH1D("hpichth","#pi^{#pm} ;#theta_{LAB}, ^#circ",1e2,0,180);
  tree_final_state->Project("hpichth","PidChargedCand.GetMomentum().Theta()*TMath::RadToDeg()");

  TH1 *hgammaph = new TH1D("hgammaph","#gamma ;#phi_{LAB}, ^#circ",1e2,-180,180);
  tree_final_state->Project("hgammaph","PidNeutralCand.GetMomentum().Phi()*TMath::RadToDeg()");
  TH1 *hpichph = new TH1D("hpichph","#pi^{#pm} ;#phi_{LAB}, ^#circ",1e2,-180,180);
  tree_final_state->Project("hpichph","PidChargedCand.GetMomentum().Phi()*TMath::RadToDeg()");

TH1 *hgammap = new TH1D("hgammap","#gamma ;|p|, GeV/c",1e2,0,10);
 tree_final_state->Project("hgammap","PidNeutralCand.GetMomentum().Mag()");
  TH1 *hpichp = new TH1D("hpichp","#pi^{#pm} ;|p|, GeV/c",1e2,0,10);
  tree_final_state->Project("hpichp","PidChargedCand.GetMomentum().Mag()");
  ////MC infor
TH1 *hmcgammath = new TH1D("hmcgammath","#gamma ;#theta_{LAB}, ^#circ",1e2,0,180);
 tree_final_state->Project("hmcgammath","MCTrack.GetMomentum().Theta()*TMath::RadToDeg()","MCTrack.fPdgCode==22");
  TH1 *hmcpichth = new TH1D("hmcpichth","#pi^{#pm} ;#theta_{LAB}, ^#circ",1e2,0,180);
  tree_final_state->Project("hmcpichth","MCTrack.GetMomentum().Theta()*TMath::RadToDeg()","abs(MCTrack.fPdgCode)==211");

  TH1 *hmcgammaph = new TH1D("hmcgammaph","#gamma ;#phi_{LAB}, ^#circ",1e2,-180,180);
  tree_final_state->Project("hmcgammaph","MCTrack.GetMomentum().Phi()*TMath::RadToDeg()","MCTrack.fPdgCode==22");
  TH1 *hmcpichph = new TH1D("hmcpichph","#pi^{#pm} ;#phi_{LAB}, ^#circ",1e2,-180,180);
  tree_final_state->Project("hmcpichph","MCTrack.GetMomentum().Phi()*TMath::RadToDeg()","abs(MCTrack.fPdgCode)==211");

TH1 *hmcgammap = new TH1D("hmcgammap","#gamma ;|p|, GeV/c",1e2,0,10);
 tree_final_state->Project("hmcgammap","MCTrack.GetMomentum().Mag()","MCTrack.fPdgCode==22");
  TH1 *hmcpichp = new TH1D("hmcpichp","#pi^{#pm} ;|p|, GeV/c",1e2,0,10);
  tree_final_state->Project("hmcpichp","MCTrack.GetMomentum().Mag()","abs(MCTrack.fPdgCode)==211");


  TCanvas c0;
  c0.Divide(3,2);
  c0.cd(1);
  hmcgammath->SetFillColor(42);
  hmcgammath->SetLineColor(1);
  hmcgammath->SetMinimum(0);
  hmcgammath->Draw();
  hgammath->SetFillColor(30);
  hgammath->SetLineColor(1);
  hgammath->SetMinimum(0);
  hgammath->Draw("same");
  c0.cd(2);
  hmcgammaph->SetFillColor(42);
  hmcgammaph->SetLineColor(1);
  hmcgammaph->SetMinimum(0);
  hmcgammaph->Draw();
  hgammaph->SetFillColor(30);
  hgammaph->SetLineColor(1);
  hgammaph->Draw("same");
  c0.cd(3);
  hmcgammap->SetFillColor(42);
  hmcgammap->SetLineColor(1);
  hmcgammap->Draw();
  hgammap->SetFillColor(30);
  hgammap->SetLineColor(1);
  hgammap->Draw("same");
  c0.cd(4);
  hmcpichth->SetFillColor(42);
  hmcpichth->SetLineColor(1);
  hmcpichth->Draw();
  hpichth->SetFillColor(30);
  hpichth->SetLineColor(1);
  hpichth->Draw("same");
  c0.cd(5);
  hmcpichph->SetFillColor(42);
  hmcpichph->SetLineColor(1);
  hmcpichph->SetMinimum(0);
  hmcpichph->Draw();
  hpichph->SetFillColor(30);
  hpichph->SetLineColor(1);
  hpichph->SetMinimum(0);
  hpichph->Draw("same");
  c0.cd(6);
  hmcpichp->SetFillColor(42);
  hmcpichp->SetLineColor(1);
  hmcpichp->Draw();
  hpichp->SetFillColor(30);
  hpichp->SetLineColor(1);
  hpichp->Draw("same");

  c0.Print("plotRes.pdf(");

  TFile *finSIG = new TFile(file_sig.Data(),"READ");
  TTree *tree_sig_raw = (TTree*)finSIG->Get("ntp3");
  TTree *tree_sig_4Cfit = (TTree*)finSIG->Get("ntp4");
  TTree *tree_sig_rho = (TTree*)finSIG->Get("ntp2");
  TTree *tree_sig_pi0 = (TTree*)finSIG->Get("ntp1");
  TTree *tree_sig_eta = (TTree*)finSIG->Get("ntp5");

  TH1 *hhcmraw = new TH1D("hhcmraw","Signal ;m, GeV",2e2,3.4,3.7);
  tree_sig_raw->Project("hhcmraw","hcm");
  TH2 *hDalitz = new TH2D("hDalitz","; m^{2}_{#pi^{+}#pi^{-}#pi^{0}}, (GeV/c^{2})^{2};m^{2}_{#pi^{+}#pi^{-}}, (GeV/c^{2})^{2}",1e3,0,10,1e3,0,10);
  tree_sig_raw->Project("hDalitz","hcd0m*hcd0m:hcd1m*hcd1m");
  TH1 *hhcmfit = new TH1D("hhcmfit","Signal ;m, GeV",1e3,2.5,4.5);
  tree_sig_4Cfit->Project("hhcmfit","hcfitm");

  TH1 *hrhom = new TH1D("hrhom","Signal ; m, GeV",1e2,0.,5.0);
  tree_sig_rho->Project("hrhom","pichpairm");
TH2 *hPeyrou_rho = new TH2D("hPeyrou_rho","Peyrou Diagram: #rho (SIG); p_{z}, GeV/c ; p_{#perp} , GeV/c",1e2,-1,7,1e3,0,2);
  tree_sig_rho->Project("hPeyrou_rho","pichpairpt:pichpairpz");
  TH1 *hpi0m = new TH1D("hpi0m","Signal ; m, GeV",1e2,0,0.3);
  tree_sig_pi0->Project("hpi0m","pi0m");
TH2 *hPeyrou_pi0 = new TH2D("hPeyrou_pi0","Peyrou Diagram: #pi^{0} (SIG); p_{z}, GeV/c; p_{#perp} , GeV/c",1e2,-1,7,1e3,0,2);
  tree_sig_pi0->Project("hPeyrou_pi0","pi0pt:pi0pz");
  TH2 *hPeyrou_hc = new TH2D("hPeyrou_hc","Peyrou Diagram: h_{c} (SIG); p_{z}, GeV/c; p_{#perp} , GeV/c",1e2,5,7,1e3,0,0.2);
  tree_sig_raw->Project("hPeyrou_hc","hcpt:hcpz");
  TH2 *hPeyrou_hc_fit = new TH2D("hPeyrou_hc_fit","Peyrou Diagram: h_{c} (SIG, 4C fit); p_{z}, GeV/c; p_{#perp} , GeV/c",1e2,5,7,1e3,0,4e-11);
  tree_sig_4Cfit->Project("hPeyrou_hc_fit","hcfitpt:hcfitpz");

TH1 *hetam = new TH1D("hetam","Signal ; m, GeV",1e2,0,5.0);
 tree_sig_eta->Project("hetam","etam");

 

  TFile *finBKG = new TFile(file_bkg.Data(),"READ");
  TTree *tree_bkg_raw = (TTree*)finBKG->Get("ntp3");
  TTree *tree_bkg_4Cfit = (TTree*)finBKG->Get("ntp4");
  TTree *tree_bkg_rho = (TTree*)finBKG->Get("ntp2");
  TTree *tree_bkg_pi0 = (TTree*)finBKG->Get("ntp1");
  TTree *tree_bkg_eta = (TTree*)finBKG->Get("ntp5");
  TH1 *hhcmraw_bkg = new TH1D("hhcmraw_bkg","Background ;m, GeV",2e2,3.4,3.7);
  tree_bkg_raw->Project("hhcmraw_bkg","hcm");
  TH2 *hDalitz_bkg = new TH2D("hDalitz_bkg","; m^{2}_{#pi^{+}#pi^{-}#pi^{0}}, (GeV/c^{2})^{2};m^{2}_{#pi^{+}#pi^{-}}, (GeV/c^{2})^{2}",1e3,0,10,1e3,0,10);
  tree_bkg_raw->Project("hDalitz_bkg","hcd0m*hcd0m:hcd1m*hcd1m");
  TH1 *hhcmfit_bkg = new TH1D("hhcmfit_bkg","Background ;m, GeV",1e3,2.5,4.5);
  tree_bkg_4Cfit->Project("hhcmfit_bkg","hcfitm");
  TH1 *hrhom_bkg = new TH1D("hrhom_bkg","Background ; m, GeV",1e2,0,5.0);
  tree_bkg_rho->Project("hrhom_bkg","pichpairm");
  TH1 *hpi0m_bkg = new TH1D("hpi0m_bkg","Background ; m, GeV",1e2,0,0.3);
  tree_bkg_pi0->Project("hpi0m_bkg","pi0m");
  TH2 *hPeyrou_pi0_bkg = new TH2D("hPeyrou_pi0_bkg","Peyrou Diagram: #pi^{0} (BKG);  p_{z}, GeV/c; p_{#perp} , GeV/c",1e2,-1,7,1e3,0,2);
  tree_bkg_pi0->Project("hPeyrou_pi0_bkg","pi0pt:pi0pz");
TH2 *hPeyrou_rho_bkg = new TH2D("hPeyrou_rho_bkg","Peyrou Diagram: #rho (BKG);  p_{z}, GeV/c; p_{#perp} , GeV/c ",1e2,-1,7,1e3,0,2);
  tree_bkg_rho->Project("hPeyrou_rho_bkg","pichpairpt:pichpairpz");
TH2 *hPeyrou_hc_bkg = new TH2D("hPeyrou_hc_bkg","Peyrou Diagram: h_{c} (SIG); p_{z}, GeV/c ; p_{#perp} , GeV/c",1e2,5,7,1e3,0,0.2);
  tree_bkg_raw->Project("hPeyrou_hc_bkg","hcpt:hcpz");
  TH2 *hPeyrou_hc_fit_bkg = new TH2D("hPeyrou_hc_fit_bkg","Peyrou Diagram: h_{c} (SIG, 4C fit); p_{z}, GeV/c; p_{#perp} , GeV/c",1e2,5,7,1e3,0,4e-11);
  tree_bkg_4Cfit->Project("hPeyrou_hc_fit_bkg","hcfitpt:hcfitpz");
TH1 *hetam_bkg = new TH1D("hetam_bkg","Signal ; m, GeV",1e2,0,5.0);
 tree_bkg_eta->Project("hetam_bkg","etam");

  TCanvas c1;
  c1.SetLogy(1);
  c1.Divide(5,2);
  c1.cd(1);
  hpi0m->Draw();
  c1.cd(2);
  hetam->Draw();
 c1.cd(3);
  hrhom->Draw();
  c1.cd(4);
  hhcmraw->Draw();
 c1.cd(5);
  hhcmfit->Draw();
  
 
  c1.cd(6);
  hpi0m_bkg->Draw();
  c1.cd(7);
  hetam_bkg->Draw();
  c1.cd(8);
  hrhom_bkg->Draw();
  c1.cd(9);
  hhcmraw_bkg->Draw();
  c1.cd(10);
  hhcmfit_bkg->Draw();
  
  

  c1.Print("plotRes.pdf(");
  TCanvas c2;
  c2.Divide(4,2);
  c2.cd(1);
  hPeyrou_pi0->Draw("colz");
  c2.cd(2);
  hPeyrou_rho->Draw("colz");
  c2.cd(3);
  hPeyrou_hc->Draw("colz");
  c2.cd(4);
  hPeyrou_hc_fit->Draw("colz");
  c2.cd(5);
  hPeyrou_pi0_bkg->Draw("colz");
  c2.cd(6);
  hPeyrou_rho_bkg->Draw("colz");
  c2.cd(7);
  hPeyrou_hc_bkg->Draw("colz");
  c2.cd(8);
  hPeyrou_hc_fit_bkg->Draw("colz");
  c2.Print("plotRes.pdf(");

  TCanvas c3;
  c3.Divide(3,2);
  c3.cd(1);
  hDalitz->Draw("colz");
  c3.cd(4);
  hDalitz_bkg->Draw("colz");
  c3.Print("plotRes.pdf)");

  TString out="plotRes.root";
  TFile *fout = new TFile(out,"RECREATE");
  hgammath->Write();
  hgammaph->Write();
  hpichth->Write();
  hpichph->Write();
  hpi0m->Write();
  hrhom->Write();
  hhcmraw->Write();
  hhcmfit->Write();
  hpi0m_bkg->Write();
  hrhom_bkg->Write();
  hhcmraw_bkg->Write();
  hhcmfit_bkg->Write();
  hPeyrou_pi0->Write();
  hPeyrou_rho->Write();
  hPeyrou_hc->Write();
  hPeyrou_hc_fit->Write();
  hPeyrou_pi0_bkg->Write();
  hPeyrou_rho_bkg->Write();
  hPeyrou_hc_bkg->Write();
  hPeyrou_hc_fit_bkg->Write();
  hDalitz_bkg->Write();
  hDalitz->Write();
  c1.Write();
  c2.Write();
  c3.Write();
  fout->Close();
}