{ 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"); }