makeStaticData()->SetParticleMass("dilepton",0.06);

{
  TH1F * histo2 = new TH1F ("histo2","IM e+ e-",80,0.58,0.146);
  makeStaticData()->SetParticleMass("dilepton",0.06);
  
  double ebeam_min = 1.1725;
  double ebeam_max = 2.9;
  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);
  
  PReaction *rr=new PReaction("_P1 = 2.2","g","p","p pi0 [dilepton [e+ e-] g]","ppi0");
  rr->Do(histo2," _x = ([dilepton])->M() ");
  rr->Print();
  rr->loop(125000);
  histo2->Draw();
}

makeStaticData()->SetParticleMass("dilepton",0.06);

#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 "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.42/src/PParticle.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.42/src/PReaction.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.42/src/PBeamSmearing.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.42/src/PAnyDistribution.h"
#include "/w/hallb/clasg12/mkunkel/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("PI0_Hists_SAID.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 *example0 = (TH1D*)fhist->Get("hist0");
  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");
  TH1D *example65 = (TH1D*)fhist->Get("hist65");
  TH1D *example66 = (TH1D*)fhist->Get("hist66");
  TH1D *example67 = (TH1D*)fhist->Get("hist67");
  TH1D *example68 = (TH1D*)fhist->Get("hist68");
  TH1D *example69 = (TH1D*)fhist->Get("hist69");
  
  //##################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("pi0,daughter,primary");
  
  //Define the range of the c.m. sampling: G12 c.m. beam energy 1.71565, 2.55105 is 1.1 - 3.0 GeV in lab
  double q_min = 1.71565;//1.71565; //;
  double q_max = 2.55105;// 2.7//;
  model->SetRange(q_min,q_max);
  //These Histograms are for the model
  
  
  //These Histograms are for the model
  //model->AddHistogram(example0,"if (_y > 1.69534 && _y <= 1.70912) value = Eval(_x); _f =value");
  model->AddHistogram(example1,"if (_y > 1.70912 && _y <= 1.72279) value = Eval(_x); _f =value");
  model->AddHistogram(example2,"if (_y > 1.72279 && _y <= 1.73635) value = Eval(_x); _f =value");
  model->AddHistogram(example3,"if (_y > 1.73635 && _y <= 1.74981) value = Eval(_x); _f =value");
  model->AddHistogram(example4,"if (_y > 1.74981 && _y <= 1.76316) value = Eval(_x); _f =value");
  model->AddHistogram(example5,"if (_y > 1.76316 && _y <= 1.77642) value = Eval(_x); _f =value");
  model->AddHistogram(example6,"if (_y > 1.77642 && _y <= 1.78957) value = Eval(_x); _f =value");
  model->AddHistogram(example7,"if (_y > 1.78957 && _y <= 1.80263) value = Eval(_x); _f =value");
  model->AddHistogram(example8,"if (_y > 1.80263 && _y <= 1.8156) value = Eval(_x); _f =value");
  model->AddHistogram(example9,"if (_y > 1.8156 && _y <= 1.82847) value = Eval(_x); _f =value");
  model->AddHistogram(example10,"if (_y > 1.82847 && _y <= 1.84126) value = Eval(_x); _f =value");
  model->AddHistogram(example11,"if (_y > 1.84126 && _y <= 1.85395) value = Eval(_x); _f =value");
  model->AddHistogram(example12,"if (_y > 1.85395 && _y <= 1.86656) value = Eval(_x); _f =value");
  model->AddHistogram(example13,"if (_y > 1.86656 && _y <= 1.87909) value = Eval(_x); _f =value");
  model->AddHistogram(example14,"if (_y > 1.87909 && _y <= 1.89153) value = Eval(_x); _f =value");
  
  model->AddHistogram(example15,"if (_y > 1.89153 && _y <= 1.90389) value = Eval(_x); _f =value");
  model->AddHistogram(example16,"if (_y > 1.90389 && _y <= 1.91617) value = Eval(_x); _f =value");
  model->AddHistogram(example17,"if (_y > 1.91617 && _y <= 1.92837) value = Eval(_x); _f =value");
  model->AddHistogram(example18,"if (_y > 1.92837 && _y <= 1.9405) value = Eval(_x); _f =value");
  model->AddHistogram(example19,"if (_y > 1.9405 && _y <= 1.95255) value = Eval(_x); _f =value");
  model->AddHistogram(example20,"if (_y > 1.95255 && _y <= 1.96453) value = Eval(_x); _f =value");
  model->AddHistogram(example21,"if (_y > 1.96453 && _y <= 1.97643) value = Eval(_x); _f =value");
  model->AddHistogram(example22,"if (_y > 1.97643 && _y <= 1.98826) value = Eval(_x); _f =value");
  model->AddHistogram(example23,"if (_y > 1.98826 && _y <= 2.00003) value = Eval(_x); _f =value");
  model->AddHistogram(example24,"if (_y > 2.00003 && _y <= 2.01172) value = Eval(_x); _f =value");
  model->AddHistogram(example25,"if (_y > 2.01172 && _y <= 2.02335) value = Eval(_x); _f =value");
  model->AddHistogram(example26,"if (_y > 2.02335 && _y <= 2.03491) value = Eval(_x); _f =value");
  model->AddHistogram(example27,"if (_y > 2.03491 && _y <= 2.0464) value = Eval(_x); _f =value");
  model->AddHistogram(example28,"if (_y > 2.0464 && _y <= 2.05783) value = Eval(_x); _f =value");
  model->AddHistogram(example29,"if (_y > 2.05783 && _y <= 2.0692) value = Eval(_x); _f =value");
  model->AddHistogram(example30,"if (_y > 2.0692 && _y <= 2.08051) value = Eval(_x); _f =value");
  model->AddHistogram(example31,"if (_y > 2.08051 && _y <= 2.09175) value = Eval(_x); _f =value");
  model->AddHistogram(example32,"if (_y > 2.09175 && _y <= 2.10293) value = Eval(_x); _f =value");
  model->AddHistogram(example33,"if (_y > 2.10293 && _y <= 2.11406) value = Eval(_x); _f =value");
  model->AddHistogram(example34,"if (_y > 2.11406 && _y <= 2.12513) value = Eval(_x); _f =value");
  model->AddHistogram(example35,"if (_y > 2.12513 && _y <= 2.13613) value = Eval(_x); _f =value");
  model->AddHistogram(example36,"if (_y > 2.13613 && _y <= 2.14709) value = Eval(_x); _f =value");
  model->AddHistogram(example37,"if (_y > 2.14709 && _y <= 2.15798) value = Eval(_x); _f =value");
  model->AddHistogram(example38,"if (_y > 2.15798 && _y <= 2.16883) value = Eval(_x); _f =value");
  model->AddHistogram(example39,"if (_y > 2.16883 && _y <= 2.17962) value = Eval(_x); _f =value");
  model->AddHistogram(example40,"if (_y > 2.17962 && _y <= 2.19035) value = Eval(_x); _f =value");
  model->AddHistogram(example41,"if (_y > 2.19035 && _y <= 2.20103) value = Eval(_x); _f =value");
  model->AddHistogram(example42,"if (_y > 2.20103 && _y <= 2.21167) value = Eval(_x); _f =value");
  model->AddHistogram(example43,"if (_y > 2.21167 && _y <= 2.22225) value = Eval(_x); _f =value");
  model->AddHistogram(example44,"if (_y > 2.22225 && _y <= 2.23278) value = Eval(_x); _f =value");
  model->AddHistogram(example45,"if (_y > 2.23278 && _y <= 2.24326) value = Eval(_x); _f =value");
  model->AddHistogram(example46,"if (_y > 2.24326 && _y <= 2.25369) value = Eval(_x); _f =value");
  model->AddHistogram(example47,"if (_y > 2.25369 && _y <= 2.26407) value = Eval(_x); _f =value");
  model->AddHistogram(example48,"if (_y > 2.26407 && _y <= 2.27441) value = Eval(_x); _f =value");
  model->AddHistogram(example49,"if (_y > 2.27441 && _y <= 2.2847) value = Eval(_x); _f =value");
  model->AddHistogram(example50,"if (_y > 2.2847 && _y <= 2.29495) value = Eval(_x); _f =value");
  model->AddHistogram(example51,"if (_y > 2.29495 && _y <= 2.30514) value = Eval(_x); _f =value");
  model->AddHistogram(example52,"if (_y > 2.30514 && _y <= 2.3153) value = Eval(_x); _f =value");
  model->AddHistogram(example53,"if (_y > 2.3153 && _y <= 2.32541) value = Eval(_x); _f =value");
  model->AddHistogram(example54,"if (_y > 2.32541 && _y <= 2.33547) value = Eval(_x); _f =value");
  model->AddHistogram(example55,"if (_y > 2.33547 && _y <= 2.34549) value = Eval(_x); _f =value");
  model->AddHistogram(example56,"if (_y > 2.34549 && _y <= 2.35547) value = Eval(_x); _f =value");
  model->AddHistogram(example57,"if (_y > 2.35547 && _y <= 2.36541) value = Eval(_x); _f =value");
  model->AddHistogram(example58,"if (_y > 2.36541 && _y <= 2.37531) value = Eval(_x); _f =value");
  model->AddHistogram(example59,"if (_y > 2.37531 && _y <= 2.38516) value = Eval(_x); _f =value");
  model->AddHistogram(example60,"if (_y > 2.38516 && _y <= 2.39498) value = Eval(_x); _f =value");
  model->AddHistogram(example61,"if (_y > 2.39498 && _y <= 2.40475) value = Eval(_x); _f =value");
  model->AddHistogram(example62,"if (_y > 2.40475 && _y <= 2.41449) value = Eval(_x); _f =value");
  model->AddHistogram(example63,"if (_y > 2.41449 && _y <= 2.42418) value = Eval(_x); _f =value");
  model->AddHistogram(example64,"if (_y > 2.42418 && _y <= 2.43384) value = Eval(_x); _f =value");
  model->AddHistogram(example65,"if (_y > 2.43384 && _y <= 2.44346) value = Eval(_x); _f =value");
  model->AddHistogram(example66,"if (_y > 2.44346 && _y <= 2.45304) value = Eval(_x); _f =value");
  model->AddHistogram(example67,"if (_y > 2.45304 && _y <= 2.46258) value = Eval(_x); _f =value");
  model->AddHistogram(example68,"if (_y > 2.46258 && _y <= 2.47209) value = Eval(_x); _f =value");
  //model->AddHistogram(example69,"if (_y > 2.47209 && _y <= 2.48156) value = Eval(_x); _f =value");
  model->AddHistogram(example69,"if (_y > 2.47209 && _y <= 2.7) value = Eval(_x); _f =value");
  
  
  //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,"_f = _f * Eval(t_lab);");
  
  
  //model->SetNpy(300);
  //model->SetNpx(300);  //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);
  //TH1F * histo2 = new TH1F ("histo2","IM e+ e-",80,0.58,0.146);

  //gROOT->Reset();
  PUtils::SetSeed(input*1000); //input
  makeStaticData()->SetParticleMass("dilepton",0.06);
  PReaction my_reaction("_P1 = 2.2","g","p","p pi0 [dilepton [e+ e-] g]","pi0_dalitz_QEDFF",1,0,0,0);

  //my_reaction.Do(histo2," _x = ([dilepton])->M() ");  

  my_reaction.Print();
  
  my_reaction.Loop(events);
  //histo2->Draw();

}