//#include "loadPluto.h";
//Program to generate multiple PLUTO root file
//Author Michael C. Kunkel & John Goetz
double beam_fn(double* x, double* par)
{
    return 1.;
}

void Try_T_Rho(){   


     //First, we create the general purpose distribution
    //model:
    PAnyDistribution* decay = 
        new PAnyDistribution("t_slope",
			     "A function to add a new t-slope");
    decay->Add("q,     parent");
    decay->Add("p,     daughter");
    decay->Add("rho0,  daughter");
    //This is the cache for the undistorted data
    //It is needed because the mandelstam t is a non-uniform
    //distribution. The size and binning of the cache must
    //be chosen such, that during runtime (better: during Preheating) the statistic
    //is sufficiently filled
    TH1F * cache  = new TH1F ("cache","Rho0 t cache",400,-4.0,.0);
    

    //For the following we we have to know that all particles (but the daughters) are in the parent rest frame
    //the daughters are in their rest frame (i.e. the parent)
    //therefore we have to boost into the parent frame
    //the parent is indicated by "_parent"
    //N.B.: "t" is reserved in TFormula, do NOT use it
    decay->AddEquation(cache,"beam = _parent->GetBeam(); beam->Boost(_parent) ; t1 = (beam - [rho0])->Mag2(); _x = t1;");
    //This is the final equation. The distribution (the probability function)
    //must be stored in "_f"
    decay->AddEquation("_f = exp( 7.1*t1 );");

    //Remember, AnyDistribution is a rejection method. Therefore
    //it can happen, that parts of the phase space, where _f has a 
    //large probability, is not well populated by the generated events
    //In this case, the event loop will run forever, as Pluto tries
    //to match the shape defined by _f.
    //The following factor is the maximum enhancement factor to avoid such
    //deadlocks.
    //N.B.: It directly scales with the computing time!!!
    decay->SetMaxEnhancementFactor(2);
    //Hint: in this configuration, the sampling of 100kEvent takes ~30min

    //Add this model to the Pluto data base:
    makeDistributionManager()->Add(decay);

    //Construct the reaction, as usual:

    double ebeam_min = 1.0;
    double ebeam_max = 5.7;
    PBeamSmearing *beam_smear = new PBeamSmearing("beam_smear", "Beam smearing");

   
    TF1* beam_smear_fn = new TF1("beam_smear_fn", "exp(12.87 - 0.6108*x+0.0387*x*x)", ebeam_min, ebeam_max);    

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

    PReaction my_reaction("_T1 = 2.2","g","p","p rho0 [pi+ pi-]","t");

    TH1F * histo1 = new TH1F ("histo1","rho0 t",100,-4,0);
    TH1F * histo3 = new TH1F ("histo3","cos theta of rho0",50,-1.,1.);

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


    my_reaction.Print();   //The "Print()" statement is optional

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