{
	   
    TH1F * histo1 = new TH1F ("histo1","Beam",100,1.,6);
    TH1F * histo2 = new TH1F ("histo2","GetBeam",100,1.,6);
    TH1F * histo3 = new TH1F ("histo3","Beam Reconstructed",100,1.,6);
    TH1F * histo4 = new TH1F ("histo4","e+e-",100,0.0,0.6);

//For Beam Smearing    
    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);
//For Form Factor =1 ie. QED 
   ((PDalitzDecay * )makeDistributionManager()->GetDistribution("eta_dalitz"))->SetUseQED(1);

//Set-up Reaction
    PReaction my_reaction("_P1 = 2.2","g","p","p eta [dilepton [e+ e-] g]","eta_dalitz",1,0,0,0);

//Do Histograms
    my_reaction.Do(histo1,"_x = [g,1]->E()");
    my_reaction.Do(histo2,"_x = [g+p]->GetBeam()");
    my_reaction.Do(histo3,"ene = (([p,2] + [eta]) -[p,1]) ;_x=ene->E()");
    my_reaction.Do(histo4,"mass = [e+] + [e-] ;_x=mass->M()");

//Do Reaction
    my_reaction.Print();
    my_reaction.Loop(10000);

//Draw    
    TCanvas c1; TCanvas c2; TCanvas c3; TCanvas c4;
    c1->cd();
    histo1->Draw();
    c2->cd();
    histo2->Draw();
    c3->cd();
    histo3->Draw();
    c4->cd();
    histo4->Draw();
}