#include "TGraphAsymmErrors.h"
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TChain.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TGraphAsymmErrors.h"
#include "TFile.h"
#include "/Users/Mike/Pluto/pluto_v5.42/src/PParticle.h"
#include "/Users/Mike/Pluto/pluto_v5.42/src/PReaction.h"
#include "/Users/Mike/Pluto/pluto_v5.42/src/PBeamSmearing.h"
#include "/Users/Mike/Pluto/pluto_v5.42/src/PAnyDistribution.h"
#include "/Users/Mike/Pluto/pluto_v5.42/plugins/scatter_mod/PScatterCrossSection.h"



void pi0_XSection(int input, int events) //int input, int events

{
//Open TFile and get histograms

  
  TFile *fhist = new TFile("/Volumes/Mac_Storage/Work_Codes/PIZERO/Dalitz_SIM/XSection/PI0_Hists.root");
  //TFile *f2hist = new TFile("/Volumes/Mac_Storage/Work_Data/DALITZ_SIMULATION/crosssectionData/Eta_CrossSection/diff_Xsection/TriggerHist.root");

//##################LOAD HISTOGRAMS HERE###################  
  TH1D *lab_profile = (TH1D*)fhist->Get("Blab_smear");
  TH1D *cm_profile = (TH1D*)fhist->Get("B_smear");
  //TH1D *trigger_profile = (TH1D*)f2hist->Get("trigger4_beam");
/*
  TH1D *example9 = (TH1D*)fhist->Get("hist9");
  
  TH1D *example10 = (TH1D*)fhist->Get("hist10");
  TH1D *example11 = (TH1D*)fhist->Get("hist11");
  TH1D *example12 = (TH1D*)fhist->Get("hist12");
  TH1D *example13 = (TH1D*)fhist->Get("hist13");
  TH1D *example14 = (TH1D*)fhist->Get("hist14");
  TH1D *example15 = (TH1D*)fhist->Get("hist15");
  TH1D *example16 = (TH1D*)fhist->Get("hist16");
  TH1D *example17 = (TH1D*)fhist->Get("hist17");
  TH1D *example18 = (TH1D*)fhist->Get("hist18");
  TH1D *example19 = (TH1D*)fhist->Get("hist19");
  
  TH1D *example20 = (TH1D*)fhist->Get("hist20");
  TH1D *example21 = (TH1D*)fhist->Get("hist21");
  TH1D *example22 = (TH1D*)fhist->Get("hist22");
  TH1D *example23 = (TH1D*)fhist->Get("hist23");
  TH1D *example24 = (TH1D*)fhist->Get("hist24");
  TH1D *example25 = (TH1D*)fhist->Get("hist25");
  TH1D *example26 = (TH1D*)fhist->Get("hist26");
  TH1D *example27 = (TH1D*)fhist->Get("hist27");
  TH1D *example28 = (TH1D*)fhist->Get("hist28");
  TH1D *example29 = (TH1D*)fhist->Get("hist29");
  
  TH1D *example30 = (TH1D*)fhist->Get("hist30");
  TH1D *example31 = (TH1D*)fhist->Get("hist31");
  TH1D *example32 = (TH1D*)fhist->Get("hist32");
  TH1D *example33 = (TH1D*)fhist->Get("hist33");
  TH1D *example34 = (TH1D*)fhist->Get("hist34");
  TH1D *example35 = (TH1D*)fhist->Get("hist35");
  TH1D *example36 = (TH1D*)fhist->Get("hist36");
 */
  TH1D *example37 = (TH1D*)fhist->Get("hist37");
  
  /*
  TH1D *example38 = (TH1D*)fhist->Get("hist38");
  TH1D *example39 = (TH1D*)fhist->Get("hist39");
  
  TH1D *example40 = (TH1D*)fhist->Get("hist40");
  TH1D *example41 = (TH1D*)fhist->Get("hist41");
  TH1D *example42 = (TH1D*)fhist->Get("hist42");
  TH1D *example43 = (TH1D*)fhist->Get("hist43");
  TH1D *example44 = (TH1D*)fhist->Get("hist44");
  TH1D *example45 = (TH1D*)fhist->Get("hist45");
*/
//##################END HISTOGRAMS HERE###################  
  
//PLUTO TIME
  
  PScatterCrossSection *model = new PScatterCrossSection("mymodel","My cross section");
  model->Add("g,grandparent,beam");
  model->Add("p,grandparent,target");
  model->Add("q,parent");
  model->Add("p,daughter");
  model->Add("eta,daughter,primary");
  
  //Define the range of the c.m. sampling: G12 c.m. beam energy 1.739,3.3268
  double q_min =2.35114; //1.70278;
  double q_max = 2.37102;//2.52442;
  model->SetRange(q_min,q_max);
//These Histograms are for the model

  
  //These Histograms are for the model
  /*
  model->AddHistogram(example9,"if (_y >1.70278 && _y <= 1.73012) value = Eval(_x); _f =value");
  
  model->AddHistogram(example10,"if (_y >1.73012 && _y <= 1.75704) value = Eval(_x); _f =value");
  
  model->AddHistogram(example11,"if (_y >1.75704 && _y <= 1.78355) value = Eval(_x); _f =value");
  
  model->AddHistogram(example12,"if (_y >1.78355 && _y <= 1.80968) value = Eval(_x); _f =value");
  
  model->AddHistogram(example13,"if (_y >1.80968 && _y <= 1.83543) value = Eval(_x); _f =value");
  
  model->AddHistogram(example14,"if (_y >1.83543 && _y <= 1.86083) value = Eval(_x); _f =value");
  
  model->AddHistogram(example15,"if (_y >1.86083 && _y <= 1.88588) value = Eval(_x); _f =value");
  
  model->AddHistogram(example16,"if (_y >1.88588 && _y <= 1.91061) value = Eval(_x); _f =value");
  
  model->AddHistogram(example17,"if (_y >1.91061 && _y <= 1.93502) value = Eval(_x); _f =value");
  
  model->AddHistogram(example18,"if (_y >1.93502 && _y <= 1.95912) value = Eval(_x); _f =value");
  
  model->AddHistogram(example19,"if (_y >1.95912 && _y <= 1.98294) value = Eval(_x); _f =value");
  
  model->AddHistogram(example20,"if (_y >1.98294 && _y <= 2.00647) value = Eval(_x); _f =value");
  
  model->AddHistogram(example21,"if (_y >2.00647 && _y <= 2.02972) value = Eval(_x); _f =value");
  
  model->AddHistogram(example22,"if (_y >2.02972 && _y <= 2.05272) value = Eval(_x); _f =value");
  
  model->AddHistogram(example23,"if (_y >2.05272 && _y <= 2.07546) value = Eval(_x); _f =value");
  
  model->AddHistogram(example24,"if (_y >2.07546 && _y <= 2.09795) value = Eval(_x); _f =value");
  
  model->AddHistogram(example25,"if (_y >2.09795 && _y <= 2.1202) value = Eval(_x); _f =value");
  
  model->AddHistogram(example26,"if (_y >2.1202 && _y <= 2.14223) value = Eval(_x); _f =value");
  
  model->AddHistogram(example27,"if (_y >2.14223 && _y <= 2.16403) value = Eval(_x); _f =value");
  
  model->AddHistogram(example28,"if (_y >2.16403 && _y <= 2.18561) value = Eval(_x); _f =value");
  
  model->AddHistogram(example29,"if (_y >2.18561 && _y <= 2.20698) value = Eval(_x); _f =value");
  
  model->AddHistogram(example30,"if (_y >2.20698 && _y <= 2.22814) value = Eval(_x); _f =value");
  
  model->AddHistogram(example31,"if (_y >2.22814 && _y <= 2.24911) value = Eval(_x); _f =value");
  
  model->AddHistogram(example32,"if (_y >2.24911 && _y <= 2.26988) value = Eval(_x); _f =value");
  
  model->AddHistogram(example33,"if (_y >2.26988 && _y <= 2.29047) value = Eval(_x); _f =value");
  
  model->AddHistogram(example34,"if (_y >2.29047 && _y <= 2.31087) value = Eval(_x); _f =value");
  
  model->AddHistogram(example35,"if (_y >2.31087 && _y <= 2.33109) value = Eval(_x); _f =value");
  
  model->AddHistogram(example36,"if (_y >2.33109 && _y <= 2.35114) value = Eval(_x); _f =value");
  */
  model->AddHistogram(example37,"if (_y >2.35114 && _y <= 2.37102) value = Eval(_x); _f =value");
  /*
  model->AddHistogram(example38,"if (_y >2.37102 && _y <= 2.39073) value = Eval(_x); _f =value");
  
  model->AddHistogram(example39,"if (_y >2.39073 && _y <= 2.41029) value = Eval(_x); _f =value");
  
  model->AddHistogram(example40,"if (_y >2.41029 && _y <= 2.42968) value = Eval(_x); _f =value");
  
  model->AddHistogram(example41,"if (_y >2.42968 && _y <= 2.44892) value = Eval(_x); _f =value");
  
  model->AddHistogram(example42,"if (_y >2.44892 && _y <= 2.46801) value = Eval(_x); _f =value");
  
  model->AddHistogram(example43,"if (_y >2.46801 && _y <= 2.48696) value = Eval(_x); _f =value");
  
  model->AddHistogram(example44,"if (_y >2.48696 && _y <= 2.50576) value = Eval(_x); _f =value");
  
  model->AddHistogram(example45,"if (_y >2.50576 && _y <= 2.52442) value = Eval(_x); _f =value");
  */
  
  
  
  //CLAS G11 data ends at c.m. of 2.84. For CLAS G12 we need energies to ~ 3.3268 so for now example64 will be used as the model to this energy range. If this is incorrect, make the cut to disregard the data
  //These Histograms are for the beam smear. There are 2 methods, 1 method for each histogram
  //Use only one method
  //Method 1: Beam profile is in c.m. and distributed as 1/E
  //Just evaluate histogram
  //model->AddEquation("b_value = (_y);");
  //model->AddHistogram(cm_profile,"beam_value = Eval(b_value);");
  //Method 2: Beam profile is in lab frame (t_lab, as in Ingo's notation from his example macro)
  //Must convert this profile to c.m. ie. t_lab = (_y*_y - g.mass*g.mass - p.mass*p.mass)/(2*p.mass*p.mass) - g.mass;
  model->AddEquation("t_lab = (_y*_y - g.mass*g.mass - p.mass*p.mass)/(2*p.mass*p.mass) - g.mass;");
  
  //model->AddHistogram(lab_profile,"beam_value = Eval(t_lab);");
  
 // model->AddHistogram(lab_profile,"_f = _f * Eval(t_lab);");
  //model->AddHistogram(trigger_profile,"_f = _f * Eval(t_lab);");

  
  //model->SetNpy(500);
  //model->SetNpx(500);  //Use this after SetNpy analysis

  makeDistributionManager()->Add(model);

//  TCanvas *q_cos = new TCanvas("q_cos","q_cos");
//  q_cos->cd();
//  model->GetFunction()->Draw("surf2");
  
//Set up reaction as usual now
  ((PDalitzDecay * )makeDistributionManager()->GetDistribution("pi0_dalitz"))->SetUseQED(1);
  
  //gROOT->Reset();
  PUtils::SetSeed(input); //input
  
  PReaction my_reaction("_P1 = 2.2","g","p","p pi0 [dilepton [e+ e-] g]","pi0_dalitz_QEDFF",1,0,0,0);
  
  TH1F * histo1 = new TH1F ("histo1","c.m.",500,q_min,q_max);  //2.1 , 2.3
  histo1->Sumw2();
  histo1->GetXaxis()->SetTitle("c.m. energy [GeV]"); 

  TH2F * histo2 = new TH2F ("histo2","c.m. vs cos_theta",500,2.1,2.3,20,-1,1);
  
  TH1F * histo4 = new TH1F ("histo4","PLUTO generated c.m. Cos#theta pi0",500,-1,1);
  histo4->GetXaxis()->SetTitle("Cos(#theta)");
  histo4->SetLineColor(kRed);
  histo4->Sumw2();
  
  //This histogram shows the beam profile:
  TH1F * beamProfile = new TH1F ("beamProfile","Beam Spectrum",500,1.1,2.9);
  beamProfile->Sumw2();
  beamProfile->GetXaxis()->SetTitle("E_{#gamma} GeV");

  my_reaction.Do(histo1,"_x = [g+p]->M()");
  my_reaction.Do(beamProfile,"_x = ([g+p]->GetBeam())->E();");
  my_reaction.Do(histo2,"_x = [g+p]->M(); mypi0 = [pi0]; mypi0->Boost([g+p]); _y = mypi0->CosTheta()");
  my_reaction.Do(histo4,"my2pi0 = [pi0]; my2pi0->Boost([g+p]); _x = my2pi0->CosTheta()");


  my_reaction.Print();

  my_reaction.Loop(events);
  
  
//  TCanvas *beams = new TCanvas("beams","beams",800,800);
//  beams->Divide(1,2);
//  beams->cd(1);
//  lab_profile->Draw();
//  beams->cd(2);
//  cm_profile->Draw();
//  
//  TCanvas *can1 = new TCanvas("can1","can1");
//  can1->cd();
//  histo1->Draw();
//
//  TCanvas *can2 = new TCanvas("can2","can2");
//  can2->cd();
//  beamProfile->Draw("ep");

//  TF1 *fb = new TF1("fb","[0]/x",1.1,5.6);
//  beamProfile->Fit("fb","R");
//  TCanvas *can3 = new TCanvas("can3","can3");
//
//  can3->cd();
//  histo2->Draw("ep");
 
  TCanvas *canAll = new TCanvas("canAll","canAll");
  canAll->cd();
  histo4->Draw();


  

  
  //can1->Print("XSECTIO_GEN_PLOTS/PLUTO_generated_cm_energy_nosmear.jpeg");
  //can1->Print("XSECTIO_GEN_PLOTS/PLUTO_generated_cm_energy_nosmear.C");

  //can2->Print("XSECTIO_GEN_PLOTS/Beam_Profile.jpeg");
  //can2->Print("XSECTIO_GEN_PLOTS/Beam_Profile.C");

  canAll->Print("XSECTIO_GEN_PLOTSPLUTO_PI0_generated_cos_theta.png");
  
  //fhist->Close();
}