#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 "/Users/Mike/Pluto/pluto_v5.40.5/src/PParticle.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PReaction.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PBeamSmearing.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PAnyDistribution.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/plugins/scatter_mod/PScatterCrossSection.h"

void eta_XSection()

{
//Open TFile and get histograms

  TFile *fhist = new TFile("/Volumes/Mac_Storage/Work_Data/DALITZ_SIMULATION/crosssectionData/Eta_CrossSection/diff_Xsection/Ntuple_EtaHist.root");

//##################LOAD HISTOGRAMS HERE###################  
  TH1D *lab_profile = (TH1D*)fhist->Get("Blab_smear");
  TH1D *cm_profile = (TH1D*)fhist->Get("B_smear");
  
  TH1D *example1 = (TH1D*)fhist->Get("hist1");
  TH1D *example2 = (TH1D*)fhist->Get("hist2");
  TH1D *example3 = (TH1D*)fhist->Get("hist3");
  TH1D *example4 = (TH1D*)fhist->Get("hist4");
  TH1D *example5 = (TH1D*)fhist->Get("hist5");
  TH1D *example6 = (TH1D*)fhist->Get("hist6");
  TH1D *example7 = (TH1D*)fhist->Get("hist7");
  TH1D *example8 = (TH1D*)fhist->Get("hist8");
  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");
  TH1D *example46 = (TH1D*)fhist->Get("hist46");
  TH1D *example47 = (TH1D*)fhist->Get("hist47");
  TH1D *example48 = (TH1D*)fhist->Get("hist48");
  TH1D *example49 = (TH1D*)fhist->Get("hist49");
  
  TH1D *example50 = (TH1D*)fhist->Get("hist50");
  TH1D *example51 = (TH1D*)fhist->Get("hist51");
  TH1D *example52 = (TH1D*)fhist->Get("hist52");
  TH1D *example53 = (TH1D*)fhist->Get("hist53");
  TH1D *example54 = (TH1D*)fhist->Get("hist54");
  TH1D *example55 = (TH1D*)fhist->Get("hist55");
  TH1D *example56 = (TH1D*)fhist->Get("hist56");
  TH1D *example57 = (TH1D*)fhist->Get("hist57");
  TH1D *example58 = (TH1D*)fhist->Get("hist58");
  TH1D *example59 = (TH1D*)fhist->Get("hist59");
  
  TH1D *example60 = (TH1D*)fhist->Get("hist60");
  TH1D *example61 = (TH1D*)fhist->Get("hist61");
  TH1D *example62 = (TH1D*)fhist->Get("hist62");
  TH1D *example63 = (TH1D*)fhist->Get("hist63");
  TH1D *example64 = (TH1D*)fhist->Get("hist64");
//##################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 = 1.739;
  double q_max = 3.3268;
  model->SetRange(q_min,q_max); 

//These Histograms are for the model
  
  //model->AddHistogram(example1,"if (_y > 1.68 && _y <= 1.69) value = Eval(_x); _f =value");
  //model->AddHistogram(example2,"if (_y > 1.69 && _y <= 1.7) value = Eval(_x); _f =value");
  //model->AddHistogram(example3,"if (_y > 1.7 && _y <= 1.71) value = Eval(_x); _f =value");
  //model->AddHistogram(example4,"if (_y > 1.71 && _y <= 1.72) value = Eval(_x); _f =value");
  //model->AddHistogram(example5,"if (_y > 1.72 && _y <= 1.73) value = Eval(_x); _f =value");
  model->AddHistogram(example6,"if (_y > 1.73 && _y <= 1.74) value = Eval(_x); _f =value");
  model->AddHistogram(example7,"if (_y > 1.74 && _y <= 1.75) value = Eval(_x); _f =value");
  model->AddHistogram(example8,"if (_y > 1.75 && _y <= 1.76) value = Eval(_x); _f =value");
  model->AddHistogram(example9,"if (_y > 1.76 && _y <= 1.77) value = Eval(_x); _f =value");
  
  model->AddHistogram(example10,"if (_y > 1.77 && _y <= 1.78) value = Eval(_x); _f =value");
  model->AddHistogram(example11,"if (_y > 1.78 && _y <= 1.79) value = Eval(_x); _f =value");
  model->AddHistogram(example12,"if (_y > 1.79 && _y <= 1.8) value = Eval(_x); _f =value");
  model->AddHistogram(example13,"if (_y > 1.8 && _y <= 1.81) value = Eval(_x); _f =value");
  model->AddHistogram(example14,"if (_y > 1.81 && _y <= 1.82) value = Eval(_x); _f =value");
  model->AddHistogram(example15,"if (_y > 1.82 && _y <= 1.83) value = Eval(_x); _f =value");
  model->AddHistogram(example16,"if (_y > 1.83 && _y <= 1.84) value = Eval(_x); _f =value");
  model->AddHistogram(example17,"if (_y > 1.84 && _y <= 1.85) value = Eval(_x); _f =value");
  model->AddHistogram(example18,"if (_y > 1.85 && _y <= 1.86) value = Eval(_x); _f =value");
  model->AddHistogram(example19,"if (_y > 1.86 && _y <= 1.87) value = Eval(_x); _f =value");
  
  model->AddHistogram(example20,"if (_y > 1.87 && _y <= 1.88) value = Eval(_x); _f =value");
  model->AddHistogram(example21,"if (_y > 1.88 && _y <= 1.89) value = Eval(_x); _f =value");
  model->AddHistogram(example22,"if (_y > 1.89 && _y <= 1.9) value = Eval(_x); _f =value");
  model->AddHistogram(example23,"if (_y > 1.9 && _y <= 1.91) value = Eval(_x); _f =value");
  model->AddHistogram(example24,"if (_y > 1.91 && _y <= 1.92) value = Eval(_x); _f =value");
  model->AddHistogram(example25,"if (_y > 1.92 && _y <= 1.93) value = Eval(_x); _f =value");
  model->AddHistogram(example26,"if (_y > 1.93 && _y <= 1.94) value = Eval(_x); _f =value");
  model->AddHistogram(example27,"if (_y > 1.94 && _y <= 1.95) value = Eval(_x); _f =value");
  model->AddHistogram(example28,"if (_y > 1.95 && _y <= 1.97) value = Eval(_x); _f =value"); //1.96-1.97
  model->AddHistogram(example29,"if (_y > 1.97 && _y <= 1.98) value = Eval(_x); _f =value");
  
  model->AddHistogram(example30,"if (_y > 1.98 && _y <= 1.99) value = Eval(_x); _f =value");
  model->AddHistogram(example31,"if (_y > 1.99 && _y <= 2.0) value = Eval(_x); _f =value");
  model->AddHistogram(example32,"if (_y > 2.0 && _y <= 2.01) value = Eval(_x); _f =value");
  model->AddHistogram(example33,"if (_y > 2.01 && _y <= 2.02) value = Eval(_x); _f =value");
  model->AddHistogram(example34,"if (_y > 2.02 && _y <= 2.03) value = Eval(_x); _f =value");
  model->AddHistogram(example35,"if (_y > 2.03 && _y <= 2.04) value = Eval(_x); _f =value");
  model->AddHistogram(example36,"if (_y > 2.04 && _y <= 2.05) value = Eval(_x); _f =value");
  model->AddHistogram(example37,"if (_y > 2.05 && _y <= 2.06) value = Eval(_x); _f =value");
  model->AddHistogram(example38,"if (_y > 2.06 && _y <= 2.07) value = Eval(_x); _f =value");
  model->AddHistogram(example39,"if (_y > 2.07 && _y <= 2.08) value = Eval(_x); _f =value");
  
  model->AddHistogram(example40,"if (_y > 2.08 && _y <= 2.09) value = Eval(_x); _f =value");
  model->AddHistogram(example41,"if (_y > 2.09 && _y <= 2.1) value = Eval(_x); _f =value");
  model->AddHistogram(example42,"if (_y > 2.1 && _y <= 2.12) value = Eval(_x); _f =value");
  model->AddHistogram(example43,"if (_y > 2.12 && _y <= 2.14) value = Eval(_x); _f =value");
  model->AddHistogram(example44,"if (_y > 2.14 && _y <= 2.16) value = Eval(_x); _f =value");
  model->AddHistogram(example45,"if (_y > 2.16 && _y <= 2.18) value = Eval(_x); _f =value");
  model->AddHistogram(example46,"if (_y > 2.18 && _y <= 2.2) value = Eval(_x); _f =value");
  model->AddHistogram(example47,"if (_y > 2.2 && _y <= 2.22) value = Eval(_x); _f =value");
  model->AddHistogram(example48,"if (_y > 2.22 && _y <= 2.24) value = Eval(_x); _f =value");
  model->AddHistogram(example49,"if (_y > 2.24 && _y <= 2.26) value = Eval(_x); _f =value");
  
  model->AddHistogram(example50,"if (_y > 2.26 && _y <= 2.28) value = Eval(_x); _f =value");
  model->AddHistogram(example51,"if (_y > 2.28 && _y <= 2.3) value = Eval(_x); _f =value");
  model->AddHistogram(example52,"if (_y > 2.3 && _y <= 2.32) value = Eval(_x); _f =value");
  model->AddHistogram(example53,"if (_y > 2.32 && _y <= 2.34) value = Eval(_x); _f =value");
  model->AddHistogram(example54,"if (_y > 2.34 && _y <= 2.36) value = Eval(_x); _f =value");
  model->AddHistogram(example55,"if (_y > 2.36 && _y <= 2.4) value = Eval(_x); _f =value");  
  model->AddHistogram(example56,"if (_y > 2.4 && _y <= 2.44) value = Eval(_x); _f =value");  
  model->AddHistogram(example57,"if (_y > 2.44 && _y <= 2.48) value = Eval(_x); _f =value");  
  model->AddHistogram(example58,"if (_y > 2.48 && _y <= 2.52) value = Eval(_x); _f =value");  
  model->AddHistogram(example59,"if (_y > 2.52 && _y <= 2.56) value = Eval(_x); _f =value"); //3.33 
  
  
  model->AddHistogram(example60,"if (_y > 2.56 && _y <= 2.6) value = Eval(_x); _f =value");
  model->AddHistogram(example61,"if (_y > 2.6 && _y <= 2.64) value = Eval(_x); _f =value");
  model->AddHistogram(example62,"if (_y > 2.64 && _y <= 2.68) value = Eval(_x); _f =value");
  model->AddHistogram(example63,"if (_y > 2.68 && _y <= 2.73) value = Eval(_x); _f =value");
  model->AddHistogram(example64,"if (_y > 2.73) value = Eval(_x); _f =value"); // && _y <= 2.84
  //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);");

  makeDistributionManager()->Add(model);

  
//Set up reaction as usual now  
  PReaction my_reaction("_P1 = 2.113","g","p","p eta [dilepton [e+ e-] g]");
  
  TH1F * histo1 = new TH1F ("histo1","c.m.",100,1.7,3.35);  //2.1 , 2.3
  histo1->Sumw2();
  histo1->GetXaxis()->SetTitle("c.m. energy [GeV]"); 

  my_reaction.Do(histo1,"_x = [g+p]->M()");
  
  TH2F * histo2 = new TH2F ("histo2","c.m. vs cos_theta",100,2.1,2.3,20,-1,1);
  my_reaction.Do(histo2,"_x = [g+p]->M(); myeta = [eta]; myeta->Boost([g+p]); _y = myeta->CosTheta()");
  
  //This histogram shows the beam profile:
  TH1F * beamProfile = new TH1F ("beamProfile","Beam Spectrum",100,1.1,5.9);
  beamProfile->Sumw2();
  beamProfile->GetXaxis()->SetTitle("E_{#gamma} GeV"); 

  my_reaction.Do(beamProfile,"_x = ([g+p]->GetBeam())->E();");
  
  
  TH1F * histo4 = new TH1F ("histo4","PLUTO generated c.m. Cos#theta",100,-1,1);
  my_reaction.Do(histo4,"my2eta = [eta]; my2eta->Boost([g+p]); _x = my2 eta->CosTheta()");
  histo4->GetXaxis()->SetTitle("Cos(#theta)"); 
  histo4->SetLineColor(kRed);
  histo4->Sumw2();
  
  my_reaction.Print();
  my_reaction.Loop(5000000);
  
  
  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();
  histo4->Draw("ep");
 
  TCanvas *canAll = new TCanvas("canAll","canAll");
  canAll->cd();
  example47->Scale(1./1000.);
  example47->Draw();
  histo4->Draw("same");

  
//  can1->Print("PLUTO_generated_cm_energy_smear.jpeg");
//  can1->Print("PLUTO_generated_cm_energy_smear.C");
//
//  can2->Print("Beam_Profile_smear.jpeg");
//  can2->Print("Beam_Profile_smear.C");
//
//  can3->Print("PLUTO_generated_cos_theta_smear.jpeg");
  
  //fhist->Close();
}