TGenPhaseSpace bug?

ROOT
1

Hi,

I am trying to generate a Dalitz plot for B+ -> K+KsKs.

So I have this very simple macro:
{
gStyle->SetPalette(1,0);
gSystem->Load(“libPhysics”);
TLorentzVector P(0,0,0,5.279);
TLorentzVector p1,p2,p3,p12,p13,p23;

Int_t n=3;
Double_t m[3] = {0.493677,0.497672,0.497672};
Double_t m12,m13,m23;

TGenPhaseSpace S;

S.SetDecay(P,n,m,“default”);

TH2F* dalitz = new TH2F(“dalitz”,“Dalitz plot”,100,0,25,100,0,25);

Double_t w;
for(int i=0;i<500000;i++) {
w = S.Generate();
p1 = S.GetDecay(0);
p2 = S.GetDecay(1);
p3 = S.GetDecay(2);
p12=p1+p2;
p13=p1+p3;
p23=p2+p3;

m12=p12.M2();
m13=p13.M2();
m23=p23.M2();

dalitz->Fill(m12,m13);

//cout<<w<<" ; "<<S.GetWtMax() <<endl;

}
dalitz->Draw(“COLZ”);
}

This should generate a perfectly flat Dalitz plot, but instead there seems to be some accumulation on the left and right side.
The plot can be seen at:
stanford.edu/~joshmt/dalitz.gif

Does anyone know what is wrong? It looks like someone asked this question a couple years ago (root.cern.ch/root/roottalk/roottalk01/2364.html) and never got an answer.

josh

2

I will answer my own question:

My code was the following:

[quote]
Double_t w;
for(int i=0;i<500000;i++) {
w = S.Generate();
p1 = S.GetDecay(0);
p2 = S.GetDecay(1);
p3 = S.GetDecay(2);
p12=p1+p2;
p13=p1+p3;
p23=p2+p3;

m12=p12.M2();
m13=p13.M2();
m23=p23.M2();

dalitz->Fill(m12,m13);

}[/quote]

To make TGenPhaseSpace work properly, the histogram must be filled using the weight returned by Generate().
When I change dalitz->Fill(m12,m13); to dalitz->Fill(m12,m13,w); then everything works fine.

Thanks to Denis Dujmic for noticing that the PhaseSpace.C tutorial provided with root shows this correct implementation. :blush:

josh