Dalitz Plot K(1270) Decay

(c) Randy Johnson - 1998

Mathematica notebook to plot the phase space for K(1270)  -> (Subscript[K, 0]) \(TraditionalForm\`\(π\^+\)\)  \(TraditionalForm\`\(π\^-\)\).
Choose axes of Dalitz plot to be E_K and E_pi.  Projections for the other π should be the same.

First define the masses:

K2mass = 1273.;
Kmass = 498.;
Pimass = 139.;

Define my to projection histograms and the Dalitz plot table:

EK = Table[0.,{100}];
Epi = Table[0.,{100}];
Epi2 = Table[0.,{100}];
Dalitz=Table[0.,{50},{50}];

Define function to add to table:

Itable[EE_,Emin_,Emax_,Ntable_] = 
	Max[Min[Floor[(EE-Emin)/(Emax-Emin)*Ntable +1.],Ntable],1]

eq1 = p3^2 == p1^2 + p2^2 + 2 p1 p2 cost

Clear[e1,e2,e3]
eq2 = eq1/.{p1->Sqrt[e1^2-Pimass^2],p2->Sqrt[e2^2-Kmass^2],
	p3->Sqrt[e3^2-Pimass^2]}

costheta=Simplify[cost/.Solve[eq2,cost][[1]]]

Monte Carlo portion of calculation.  First select E_Pi (e1) and E_K (e2)  randomly within the limits of energy conservation.  Cacluate the energy of the other pion.  Then check if energy and momentum conservation holds.  First, check to see if the energy of the  pion is greater than a pion mass.  Then check to see if momentum can be conserved (this translates to the cosine of the angle between the first pion and the kaon is within the real limits of -1<=cos(θ)<=1.  If it is, increment the histograms.  On the Sparc 10, this procedure takes about 3 seconds per 100 events.

n=0;
e3=K2mass;
While[n<20000,
	e1 = (K2mass-Kmass-2*Pimass)*Random[]+Pimass;
	e2 = (K2mass-Kmass-2*Pimass)*Random[]+Kmass;
	e3 = K2mass -e1 -e2;
	If[e3>=Pimass,
		If[Abs[costheta]<=1,
			i=Itable[e1,0.,K2mass,100];
			j=Itable[e2,0.,K2mass,100];
			k=Itable[e3,0.,K2mass,100];
			Epi[[i]] += 1.; 
			EK[[j]] += 1.;
			Epi2[[k]] += 1.;
			i=Itable[e1,0.,K2mass,50];
			j=Itable[e2,0.,K2mass,50];
			Dalitz[[j,i]] +=1.;
			n++;
			If[Floor[n/1000]==n/1000,
				Print["Events Generated = ",n]]
		];
	];
]//Timing

Now plot the first pion energy distribution:

plotpi1=ListPlot[Table[{(i-.5)/100*K2mass,Epi[[i]]},{i,1,100}],
	PlotJoined->True]

Check to see that we get all the events:

Apply[Plus,Epi]

Plot the second pions energy distribution.  To within statistics, this should be the same as the first:

plotpi2=ListPlot[Table[{(i-.5)/100*K2mass,Epi2[[i]]},{i,1,100}],
	PlotJoined->True]

Show[plotpi1,plotpi2]

Plot the kaon's energy distribution:

ListPlot[Table[{(i-.5)/100*K2mass,EK[[i]]},{i,1,100}],
	PlotJoined->True]

Do a contour plot of the pion energy vs. the kaon energy:

plot1 =ContourPlot[Dalitz[[Itable[y,0,K2mass,50],
	Itable[x,0,K2mass,50]]]
	,{x,0,K2mass},{y,0,K2mass}
	,ContourShading->False,PlotPoints->50]

Add the energy conservation lines to the plot:

plot2=Show[plot1,
	Graphics[Line[{{Pimass,Kmass},
	{Pimass,(K2mass-2*Pimass)},
	{(K2mass-Pimass-Kmass),Kmass},
	{Pimass,Kmass}}]],
	PlotRange->{{0,K2mass},{0,K2mass}}
	]

Note:  the contour lines outside the arise from the bin resolution of our original table.