//#include "loadPluto.h";
//Program to generate multiple PLUTO root file
//Author Michael C. Kunkel
#include "TH1.h"
#include "TH2.h"
#include "TH3.h"
#include "TChain.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"

#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.40/src/PTCrossWeight.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.40/src/PChannelModel.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.40/src/PParticle.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.40/src/PReaction.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.40/src/PBeamSmearing.h"
#include "/w/hallb/clasg12/mkunkel/PLUTO/pluto_v5.40/src/PAnyDistribution.h"


void SIMULATE_Eta_Dalitz_CROSSSECTION_2(){
    
//ADDED t-slope rejector here
	PAnyDistribution* decay = new PAnyDistribution("t_slope","A function to add a new t-slope");
    decay->Add("q,     parent");
    decay->Add("p,     daughter");
    decay->Add("eta,  daughter");

    TH1F * cache  = new TH1F ("cache","eta t cache",400,-4.0,.0);
    
    decay->AddEquation(cache,"beam = _parent->GetBeam(); beam->Boost(_parent) ; t1 = (beam - [eta])->Mag2(); _x = t1;");

    decay->AddEquation("_f = exp( 2.57*t1 );");  //2.75609 was from data fit

    decay->SetMaxEnhancementFactor(10);

    makeDistributionManager()->Add(decay);
 
//Here is Cross Section information
    if(!makeStaticData()->IsParticleValid("g + p")) {
	makeStaticData()->AddParticle(14001, "g + p",
				      makeStaticData()->GetParticleMass("p"));
	makeStaticData()->AddAlias("g + p","g+p");
    }
    
    double x_gp_etaval[] = { 774.0, 824.0, 874.0, 924.0, 975.0, 1025.0, 1073.0, 1124.0, 1175.0, 1225.0, 1277.0, 1326.0, 1374.0, 1429.0, 1480.0, 1529.0, 1575.0, 1626.0, 1674.0, 1721.0, 1776.0, 1829.0, 1878.0, 1930.0, 1978.0, 2025.0, 2073.0, 2123.0, 2174.0, 2225.0, 2277.0, 2351.0, 2450.0, 2550.0, 2700.0, 2887.0 };
    double x_gp_etaerrminus[] = { 24.0, 24.0, 24.0, 24.0, 25.0, 25.0, 23.0, 24.0, 25.0, 25.0, 27.0, 26.0, 24.0, 29.0, 30.0, 29.0, 25.0, 26.0, 24.0, 21.0, 26.0, 29.0, 28.0, 30.0, 28.0, 25.0, 23.0, 23.0, 24.0, 25.0, 27.0, 51.0, 50.0, 50.0, 100.0, 87.0 };
    double x_gp_etaerrplus[] = { 26.0, 26.0, 26.0, 26.0, 25.0, 25.0, 27.0, 26.0, 25.0, 25.0, 23.0, 24.0, 26.0, 21.0, 20.0, 21.0, 25.0, 24.0, 26.0, 29.0, 24.0, 21.0, 22.0, 20.0, 22.0, 25.0, 27.0, 27.0, 26.0, 25.0, 23.0, 49.0, 50.0, 50.0, 100.0, 113.0 };
    double y_gp_etaval[] = { 16.3005, 14.8474, 10.6035, 6.94188, 4.2497, 3.07808, 3.25968, 3.166, 3.09223, 2.82842, 2.82017, 2.60072, 2.66501, 2.07136, 1.70988, 1.85793, 1.77559, 1.56821, 1.59037, 1.47269, 1.47374, 1.40882, 1.45883, 1.48256, 1.28267, 1.18375, 1.20839, 1.31939, 1.38012, 1.19128, 0.934107, 0.977936, 0.816615, 0.758357, 0.739331, 0.633444 };
    double y_gp_etaerrminus[] = { 1.05282, 0.825122, 0.608457, 0.392925, 0.243394, 0.178359, 0.183611, 0.175687, 0.172351, 0.162073, 0.164635, 0.421538, 0.436161, 0.340327, 0.280763, 0.305414, 0.29266, 0.25361, 0.278256, 0.249551, 0.250321, 0.242164, 0.257068, 0.257329, 0.226804, 0.21209, 0.238498, 0.245857, 0.275481, 0.232461, 0.234018, 0.191841, 0.144086, 0.138337, 0.130605, 0.157451 };
    double y_gp_etaerrplus[] = { 1.05282, 0.825122, 0.608457, 0.392925, 0.243394, 0.178359, 0.183611, 0.175687, 0.172351, 0.162073, 0.164635, 0.421538, 0.436161, 0.340327, 0.280763, 0.305414, 0.29266, 0.25361, 0.278256, 0.249551, 0.250321, 0.242164, 0.257068, 0.257329, 0.226804, 0.21209, 0.238498, 0.245857, 0.275481, 0.232461, 0.234018, 0.191841, 0.144086, 0.138337, 0.130605, 0.157451 };
    double y_gp_etastatminus[] = { 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, 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, 0.0, 0.0 };
    double y_gp_etastatplus[] = { 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, 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, 0.0, 0.0 };
    int gp_eta_numpoints = 36;

    
    double thr=makeStaticData()->GetParticleMass("g") + makeStaticData()->GetParticleMass("p")+ makeStaticData()->GetParticleMass("eta");
    //Better is to convert to the SECONDARY_MODELS standard

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


    //TGraph *gp_eta_gra = new TGraph(14,x_pn_eta,y_pn_eta);
    TGraphAsymmErrors *gp_eta_gra = new TGraphAsymmErrors(gp_eta_numpoints, x_gp_etaval, y_gp_etaval, x_gp_etaerrminus, x_gp_etaerrplus, y_gp_etaerrminus, y_gp_etaerrplus);
    PTCrossWeight * gp_eta_cross = new PTCrossWeight("g + p_to_p_eta/tcross", "Eta cross section in the g+p reaction",-1);
    gp_eta_cross->Add("g,grandparent,beam");
    gp_eta_cross->Add("p,grandparent,target");
    gp_eta_cross->Add("q,parent");
    gp_eta_cross->Add("p,daughter");
    gp_eta_cross->Add("eta,daughter");
    gp_eta_cross->SetCrossSection(gp_eta_gra,kTRUE); 
    gp_eta_cross->SetScaling(.000001); //convert from [ub] to [b]
    makeDistributionManager()->Add(gp_eta_cross);

//#####END CROSS SECTION INFORMATION
    
//Set up reaction as usual	

    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); //*(0.288*10^(-3) + 0.555*10^(-2)*(x^(-4.59)))

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

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

    gROOT->Reset();

    char nam1[80] = "/volatile/clas/clasg12/mkunkel/ETA_SIM/PLUTOGEN_ROOT_FILES/eta_";
    
    char nam2[15] = "_257dalitz";
   
    for(int ij = 51; ij<=100; ij++){
		PUtils::SetSeed(ij);  //by default system time is used
		char c[10];
		sprintf(c,"%d",ij);
		//string creater = nam1 + c + nam2;
		char creater[75];
		sprintf(creater,"%s%s%s",nam1,c,nam2);
		cout << creater<<endl;

    PReaction my_reaction("_P1 = 2.2","g","p","p eta [dilepton [e+ e-] g]",creater,1,0,0,0);
    
    
    //TH1F * histo1 = new TH1F ("histo1","eta t",100,-4,0);
    //TH1F * histo3 = new TH1F ("histo3","cos theta of eta",50,-1.,1.);

    //my_reaction.Do(histo1,"beam = [g+p]->GetBeam(); t1 = (beam - [eta])->Mag2(); _x = t1;");
    //my_reaction.Do(histo3,"_eta=[eta]; _eta->Boost([g+p]); _x= cos(_eta->Theta())");

    //Make a dummy loop to fill the AnyDistribution with some statistics:
    my_reaction.Preheating(100);

    //my_reaction.Output("eta_dalitz.lst","px=[e+]->Px(); py=[e+]->Py(); pz=[e+]->Pz();\n; echo $px,$py,$pz");    
    //_T1 is the beam kinetic energy
    //one can use _P1 for the beam momentum
    //and _T2, _P2 in addition for collider experiments (separated by ';')

    //This histogram shows the beam profile:
    //TH1F * beamProfile = new TH1F ("beamProfile","Pz of beam",100,0,6);
    //my_reaction.Do(beamProfile,"_x = [g + p]->Pz();");

    // my_reaction.Print();   //The "Print()" statement is optional
    my_reaction.Loop(250000);
    
    //beamProfile->Draw();
    /*
    TCanvas *c1 = new TCanvas("c1","c1");
    c1->cd();   
    histo1->Draw();
    TCanvas *c2 = new TCanvas("c2","c2");
    c2->cd();
    beamProfile->Draw();
    */
    //my_reaction.Close(); 
    }
    
}