{
    PData::SetWeightVersion(0); //Disable 1/N_ev * BR weighting
    //First define our histograms
    TH2F *histo  = new TH2F ("histo",  "Dalitz", 30, 1., 4., 30, 1., 4.);
    TH1F *histo1 = new TH1F ("histo1", "p/pi0 missing mass", 100, 0.1, 2.0);
    TH1F *histo2 = new TH1F ("histo2", "cos theta of pp", 20, -1., 1.);
    TH1F *histo3 = new TH1F ("histo3", "cos theta of D+ decay", 20, -1., 1.);
    TH2F *histo_p1 = new TH2F ("histo_p1", "px vs. vy of the pp pair", 30, -1., 1., 30, -1., 1.);
    //If you want another decay angle of the D, uncomment this line:
    ((PScatterDistribution *)makeDistributionManager()->GetDistribution("pp_delta_waves1"))->SetAngleFunction(new TF1("delme","(1+3*x*x)*.25",-1,1));
 
    //No D+ angular distribution?
    //makeDistributionManager()->Disable("pp_delta+_angle");
    //Define the reaction as usual
    PReaction my_reaction(3.13, "p", "p", "p D+ [p pi0]", "delme", 1, 0, 0, 0);
    //Combine the masses of p,1 and pi0, p,2, pi0 to the Dalitz plot
    my_reaction.Do(histo, "_x = ([p,1] + [pi0])->M2() ; _y = ([p,2] + [pi0])->M2() ");
    //Missing mass of the p2 and pi0 pair:
    my_reaction.Do(histo1, "miss= [p + p]- ( [p,2]+ [pi0] );_x=miss->M()");
    //pp aligment angle: Important: copy the proton before boosting it, otherwise you
    //will have the boosted particle on tape
    my_reaction.Do(histo2, "p1=[p,1]; p1->Boost([p + p]); _x=cos(p1->Theta())");
    //The decay angle of the D+:
    my_reaction.Do(histo3, "_pi0=[pi0]; _d=[D+]; _pi0->Rot(_d); _d->Rot(_d); _pi0->Boost(_d); _x= cos(_pi0->Theta())");
    //Momentum of a pair
    //All methods of a PParticle of the type ->XXX(void) can be used
    //Syntax can be ->XXX() or .XXX()  (as you like)
    my_reaction.Do(histo_p1, "pp_pair= [p,1] + [p,2];_x=pp_pair.Px(); _y=pp_pair->Py()");
    my_reaction.Loop(50000);
    TCanvas *c1 = new TCanvas();
    c1->Divide(2,2);
    c1->cd(1); 
    gStyle->SetPalette(1,0);
    histo->Draw("colz");
    c1->cd(2); 
    histo2->Draw("e1");
    c1->cd(3); 
    histo1->Draw("e1");
    c1->cd(4); 
    histo3->Draw("e1");
    TCanvas *c2 = new TCanvas();
    histo_p1->Draw("colz");
}