model->AddHistogram(example,"value = Eval(_x); _f = _y * value");

model->SetRange(2.2,2.22); //in GeV

my_reaction.Do(histo1,"_x = [g+p]->M()");

my_reaction.Do(beamProfile,"_x = ([g+p]->GetBeam())->E();");

#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PAngularDistribution.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PProjector.h"
#include "/Users/Mike/Pluto/pluto_v5.40.5/src/PF2EvalBatch.h"

#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"


// 2.2 - 2.22 

void eta_XSection()
{
	
  
  //From data on eta differential cross section as function of c.m. cos(theta) for c.m. energy 2.2 -> 2.22 GeV
  //Making a TGraph then fit the TGraph to a spline
  const int p7694_d47x1y1_numpoints = 17;

  double p7694_d47x1y1_xval[p7694_d47x1y1_numpoints] = { -0.84, -0.75, -0.65, -0.55, -0.45, -0.35, -0.25, -0.15, -0.05, 0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75 };
	double p7694_d47x1y1_xerrminus[p7694_d47x1y1_numpoints] = { 0.040000000000000036, 0.050000000000000044, 0.04999999999999993, 0.04999999999999993, 0.04999999999999999, 0.050000000000000044, 0.04999999999999999, 0.05000000000000002, 0.05, 0.05, 0.04999999999999999, 0.04999999999999999, 0.04999999999999999, 0.04999999999999999, 0.050000000000000044, 0.050000000000000044, 0.050000000000000044 };
	double p7694_d47x1y1_xerrplus[p7694_d47x1y1_numpoints] = { 0.039999999999999925, 0.050000000000000044, 0.050000000000000044, 0.050000000000000044, 0.04999999999999999, 0.04999999999999999, 0.04999999999999999, 0.04999999999999999, 0.05, 0.05, 0.05000000000000002, 0.04999999999999999, 0.050000000000000044, 0.04999999999999999, 0.04999999999999993, 0.04999999999999993, 0.050000000000000044 };
	double p7694_d47x1y1_yval[p7694_d47x1y1_numpoints] = { 0.0324, 0.0374, 0.0276, 0.0233, 0.0173, 0.0122, 0.0122, 0.0169, 0.0277, 0.0333, 0.0359, 0.0375, 0.0439, 0.0625, 0.093, 0.132, 0.1634 };
	double p7694_d47x1y1_yerrminus[p7694_d47x1y1_numpoints] = { 0.005, 0.00291, 0.00185, 0.00166, 0.00139, 0.00102, 9.7E-4, 0.00117, 0.00178, 0.00205, 0.00215, 0.00217, 0.00258, 0.00386, 0.005, 0.0069, 0.0093 };
	double p7694_d47x1y1_yerrplus[p7694_d47x1y1_numpoints] = { 0.005, 0.00291, 0.00185, 0.00166, 0.00139, 0.00102, 9.7E-4, 0.00117, 0.00178, 0.00205, 0.00215, 0.00217, 0.00258, 0.00386, 0.005, 0.0069, 0.0093 };
	double p7694_d47x1y1_ystatminus[p7694_d47x1y1_numpoints] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
	double p7694_d47x1y1_ystatplus[p7694_d47x1y1_numpoints] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
	TGraphAsymmErrors *p7694_d47x1y1 = new TGraphAsymmErrors(p7694_d47x1y1_numpoints, p7694_d47x1y1_xval, p7694_d47x1y1_yval, p7694_d47x1y1_xerrminus, p7694_d47x1y1_xerrplus, p7694_d47x1y1_yerrminus, p7694_d47x1y1_yerrplus);
	p7694_d47x1y1->SetTitle("Differential Cross Section for #eta c.m. energy 2.2 GeV#rightarrow 2.22 GeV"); ///HepData/7694/d47x1y1
  p7694_d47x1y1->GetXaxis()->SetTitle("Cos(#theta)"); 
	p7694_d47x1y1->GetYaxis()->SetTitle("#frac{d#sigma}{d#Omega}"); 
  TCanvas *TgCan = new TCanvas("TgCan","TgCan");
  TgCan->cd();
	p7694_d47x1y1->Draw("AP");
  
  TSpline5 *s3 = new TSpline5("s3",p7694_d47x1y1);
  s3->Draw("same");

  
  //Now from the spline, create a TH1D histogram for PLUTO to model the Xsection with
  TH1D *example = new TH1D("example","try this",2001,-1,1);
  example->SetTitle("Extrapolated Differential Cross Section for #eta c.m. energy 2.2 GeV#rightarrow 2.22 GeV"); ///HepData/7694/d47x1y1
  example->GetXaxis()->SetTitle("Cos(#theta)"); 
	example->GetYaxis()->SetTitle("#frac{d#sigma}{d#Omega}"); 
  example->Sumw2();
  int hold = 1; //place value for histogram bin filling
//############Filling Histogram#########################  
  for (double j = -1.; j<=1.; j = j +0.001) {
    double bin_content = s3->Eval(double(j));
    if(bin_content < 0.0){bin_content = 0.0;}
    //cout<<bin_content<<endl;
    example->SetBinContent(hold,bin_content);
    hold++;
  }
//################################################  
  TCanvas *canf = new TCanvas("canf","canf");
  canf->cd();
  example->Draw("ep");
  
  
//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:
  model->SetRange(2.2,2.22);  //CLAS G12 c.m. model->SetRange(1.739,3.3268);
  

  //model->AddEquation("my_x = s3->Eval(_x);;_f = Eval(my_x)");
  model->AddHistogram(example,"value = Eval(_x); _f =_y * value");
  makeDistributionManager()->Add(model);
  
  
  
  
  PReaction my_reaction("_P1 = 2.113","g","p","p eta [dilepton [e+ e-] g]");
  
  TH1F * histo1 = new TH1F ("histo1","c.m.",100,2.1,2.5);
  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,1.0,3.0,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,2.0,2.7);
  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,"myeta->Boost([g+p]); _x = myeta->CosTheta()");
  histo4->GetXaxis()->SetTitle("Cos(#theta)"); 
  histo4->Sumw2();
  
  
  my_reaction.Print();
  
  
  my_reaction.Loop(100000);
  
  
  TCanvas *can1 = new TCanvas("can1","can1");
  can1->cd();
  histo1->Draw();
  
  TCanvas *can2 = new TCanvas("can2","can2");
  can2->cd();
  beamProfile->Draw("ep");
  TCanvas *can3 = new TCanvas("can3","can3");

  can3->cd();
  histo4->Draw("ep");
  
  //TgCan->Print("XSection_Data.jpeg");
  //canf->Print("XSection_Interpolated.jpeg");
  //can1->Print("PLUTO_generated_cm_energy.jpeg");
  //can2->Print("Beam_Profile.jpeg");
  //can3->Print("PLUTO_generated_cos_theta.jpeg");
}