double ebeam_min = 1.1725;
    double ebeam_max = 5.44575;
    PBeamSmearing *beam_smear = new PBeamSmearing("beam_smear", "Beam smearing");

   
    
    TF1* beam_smear_fn = new TF1("beam_smear_fn", "1./x", ebeam_min, ebeam_max);

    beam_smear->SetReaction("g + p");
    beam_smear->SetMomentumFunction(beam_smear_fn);
    makeDistributionManager()->Add(beam_smear);

    ((PDalitzDecay * )makeDistributionManager()->GetDistribution("eta_dalitz"))->SetUseQED(1);


    PReaction my_reaction("_P1 = 2.2","g","p","p eta [dilepton [e+ e-] g]",creater,1,0,0,0);

B_min    B_Max   W_min  W_max    x_min   x_max   x_mean  y_min  y_max    y_mean
1.03491	1.05287	1.68	1.69	-0.85	-0.9	-0.8	0.1899	0.1464	0.2334
1.03491	1.05287	1.68	1.69	-0.75	-0.8	-0.7	0.2513	0.2244	0.2782
1.03491	1.05287	1.68	1.69	-0.65	-0.7	-0.6	0.2672	0.2486	0.2858
1.03491	1.05287	1.68	1.69	-0.55	-0.6	-0.5	0.2554	0.2407	0.2701
1.03491	1.05287	1.68	1.69	-0.45	-0.5	-0.4	0.2654	0.2511	0.2797
1.03491	1.05287	1.68	1.69	-0.35	-0.4	-0.3	0.2587	0.2465	0.2709
1.03491	1.05287	1.68	1.69	-0.25	-0.3	-0.2	0.29	0.2769	0.3031
1.03491	1.05287	1.68	1.69	-0.15	-0.2	-0.1	0.284	0.2712	0.2968
1.03491	1.05287	1.68	1.69	-0.05	-0.1	0	0.2828	0.2699	0.2957
1.03491	1.05287	1.68	1.69	0.05	0	0.1	0.2699	0.2572	0.2826
1.03491	1.05287	1.68	1.69	0.15	0.1	0.2	0.3017	0.2867	0.3167
1.03491	1.05287	1.68	1.69	0.25	0.2	0.3	0.2812	0.264	0.2984
1.03491	1.05287	1.68	1.69	0.35	0.3	0.4	0.236	0.2233	0.2487
1.03491	1.05287	1.68	1.69	0.45	0.4	0.5	0.2449	0.2297	0.2601
1.03491	1.05287	1.68	1.69	0.55	0.5	0.6	0.2633	0.2455	0.2811
1.03491	1.05287	1.68	1.69	0.65	0.6	0.7	0.2439	0.2221	0.2657

for (int i=0;i<36;i++) {
	x_gp_etaval[i]*=0.001; //in GeV 
	x_gp_etaval[i]+=thr;
    }

0.774 <--Before    After-->       2.25972
0.824 <--Before    After-->       2.30972
0.874 <--Before    After-->       2.35972
0.924 <--Before    After-->       2.40972
0.975 <--Before    After-->       2.46072
1.025 <--Before    After-->       2.51072
1.073 <--Before    After-->       2.55872
1.124 <--Before    After-->       2.60972
1.175 <--Before    After-->       2.66072
1.225 <--Before    After-->       2.71072
1.277 <--Before    After-->       2.76272
1.326 <--Before    After-->       2.81172
1.374 <--Before    After-->       2.85972
1.429 <--Before    After-->       2.91472
1.48 <--Before    After-->       2.96572
1.529 <--Before    After-->       3.01472
1.575 <--Before    After-->       3.06072
1.626 <--Before    After-->       3.11172
1.674 <--Before    After-->       3.15972
1.721 <--Before    After-->       3.20672
1.776 <--Before    After-->       3.26172
1.829 <--Before    After-->       3.31472
1.878 <--Before    After-->       3.36372
1.93 <--Before    After-->       3.41572
1.978 <--Before    After-->       3.46372
2.025 <--Before    After-->       3.51072
2.073 <--Before    After-->       3.55872
2.123 <--Before    After-->       3.60872
2.174 <--Before    After-->       3.65972
2.225 <--Before    After-->       3.71072
2.277 <--Before    After-->       3.76272
2.351 <--Before    After-->       3.83672
2.45 <--Before    After-->       3.93572
2.55 <--Before    After-->       4.03572
2.7 <--Before    After-->       4.18572
2.887 <--Before    After-->       4.37272

	double xval[] = { -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 yval[] = { 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 };

    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");
}

Error: cannot open file "PF2EvalBatch.h"  /Users/Mike/Pluto/pluto_v5.40.5/plugins/scatter_mod/PScatterCrossSection.h:15:

if (npy>0) pf2->SetNpx(npy);

    //Now add the histogram to the model class, and define an equation
    //Input:  _x is cos(theta), _y is the c.m. energy
    //Output: _f: cross section

Input : _x s cos(theta), _y is differential cross section
      Output : cross section

model->AddHistogram(distribution1,"_f = Eval(_x. _y);");
model->AddHistogram(distribution2,"_f = _f * Eval(_x. _y);");

model->AddHistogram(distribution3,"_f = _f * Eval(_y);");

model->AddEquation("t_lab = (_y*_y - g.mass*g.mass - p.mass*p.mass)/(2*p.mass*p.mass) - g.mass;");
model->AddHistogram(profile,"_f = _f * Eval(t_lab);");

Input : _x s cos(theta), _y is differential cross section
      Output : _f cross section

  model1->SetRange(1.77,1.8); 
...
...
...
  model64->SetRange(2.56,2.6);

  model1->AddHistogram(example1,"value = Eval(_x); _f =_y * value");
  makeDistributionManager()->Add(model1);
  ...
  ...
  ...
  model64->AddHistogram(example64,"value = Eval(_x); _f =_y * value");
  makeDistributionManager()->Add(model64);

_x = Cos(theta)
_y = Differential Cross section

    //Input:  _x is cos(theta), _y is the c.m. energy
    //Output: _f: cross section

    TH1F *distribution = new TH1F("distribution", "Angular distribution", 10, -1 , 1 );
    distribution->SetBinContent(1,20.);
    distribution->SetBinContent(2,16.);
    distribution->SetBinContent(3,11.);
    distribution->SetBinContent(4,8.);
    distribution->SetBinContent(5,5.);
    distribution->SetBinContent(6,4.);
    distribution->SetBinContent(7,3.);
    distribution->SetBinContent(8,2.5);
    distribution->SetBinContent(9,2.);
    distribution->SetBinContent(10,1.);

    TH1F *distribution2 = new TH1F("distribution2", "Angular distribution2", 10, -1 , 1 );
    distribution2->SetBinContent(1,10.);
    distribution2->SetBinContent(2,11.);
    distribution2->SetBinContent(3,12.);
    distribution2->SetBinContent(4,13);
    distribution2->SetBinContent(5,14);
    distribution2->SetBinContent(6,17);
    distribution2->SetBinContent(7,30);
    distribution2->SetBinContent(8,40);
    distribution2->SetBinContent(9,45);
    distribution2->SetBinContent(10,60);
    
    TH1F *distribution3 = new TH1F("distribution3", "Angular distribution2", 10, -1 , 1 );
    distribution3->SetBinContent(1,1.);
    distribution3->SetBinContent(2,1.);
    distribution3->SetBinContent(3,1.);
    distribution3->SetBinContent(4,1);
    distribution3->SetBinContent(5,1);
    distribution3->SetBinContent(6,1);
    distribution3->SetBinContent(7,3);
    distribution3->SetBinContent(8,4);
    distribution3->SetBinContent(9,4);
    distribution3->SetBinContent(10,6);

    model->AddHistogram(distribution,"if (_y < 1.601) _f = Eval(_x);");
    model->AddHistogram(distribution2,"if (_y > 1.600 && _y < 1.801) _f = Eval(_x);");
    model->AddHistogram(distribution3,"if (_y > 1.800) _f = Eval(_x);");

PF2EvalBatch *GetFunction(void) {return pf2;};

model->GetFunction()->Draw("surf");

Error: Can't call PScatterCrossSection::GetFunction() in current scope Draw_q_costheta.C:206:
Possible candidates are...
(in PScatterCrossSection)
(in PAngularDistribution)
(in PDistribution)
(in TF1)
(in TFormula)
Error: non class,struct,union object GetFunction() used with . or -> Draw_q_costheta.C:206:
*** Interpreter error recovered ***

PScatterCrossSection::PScatterCrossSection(const Char_t *id, const Char_t *de) :
    PAngularDistribution(id, de) {

    beam = target = NULL;
    pf2  = NULL;
    qmax = 0.0;
    qmin = -1;
    npx = npy = -1;

    TF1 * dummy = new TF1("dummy","1",-1,1);
    SetAngleFunction(dummy);
 
} ;

model->AddEquation("t_lab = (_y*_y - g.mass*g.mass - p.mass*p.mass)/(2*p.mass*p.mass) - g.mass;");
model->AddHistogram(brem_profile,"_f = _f * Eval(t_lab);");
model->AddHistogram(trigger_profile,"_f = _f * Eval(t_lab);");