Dalitz Plot K(1270) Decay

(c) Randy Johnson - 1998

Mathematica notebook to plot the phase space for K(1270) ->K0 pi+  pi-.
Choose axes of Dalitz plot to be E_pi+ and E_pi-.  

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) Ntable)/(Emax - Emin)) + 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) - (Pimass^2)]), 
          p3 -> (Sqrt[(e3^2) - (Kmass^2)])}

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

Determine K* mass from e1:

mkstar = (Sqrt[-2.  K2mass e1 + (K2mass^2) + (Pimass^2)])

Determine weight from Breit-Wigner:

weightkstar = 
      ((gamma^2)/(4.  ((((mkstar - mm))^2) + ((gamma^2)/4)))) /. {
          gamma -> 50. , mm -> 892. }

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

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

Now plot the negative pion energy distribution -- it should look like the Breit Wigner:

plotpi1 = 
      ListPlot[Table[{(i -  .5) 5.  + 100. , 
            Epi[[i]]}, {i, 1, 100}], 
        PlotJoined -> True, PlotRange -> All,  
        PlotLabel -> "π- energy distribution"]

Plot the positive pions energy distribution:

plotpi2 = 
      ListPlot[Table[{(i -  .5) 5.  + 100. , 
            Epi2[[i]]}, {i, 1, 100}], 
        PlotJoined -> True, PlotRange -> All, 
        PlotLabel -> "π+ Energy distribution"]

Compare the two:

Show[plotpi1, plotpi2]

Plot the kaon's energy distribution:

ListPlot[
      Table[{450.  + (i -  .5) 5, 
          EK[[i]]}, {i, 1, 100}], 
      PlotJoined -> True,  PlotLabel -> "K energy distribution"]

Do a contour plot of the negative pion energy vs. positive pion energy:

plot1 = 
      ContourPlot[
        Dalitz[[Itable[y, 100. , 600. , 50], 
          Itable[x, 100. , 600. , 50]]], {x, 100. , 
          600. }, {y, 100. , 600. }, ContourShading -> False, 
        PlotPoints -> 50, PlotRange -> All, 
        AxesLabel -> {"π- energy", "π+ energy"}]

Add the energy conservation lines to the plot:

plot2 = 
      Show[plot1, 
        Graphics[
          Line[{{Pimass, Pimass}, {Pimass, K2mass - Pimass - Kmass}, {
                K2mass - Pimass - Kmass, Pimass}, {Pimass, Pimass}}]], 
        PlotRange -> {{100, 600}, {100, 600}}]

Nowcalculate kinematic boundary and plot as in problem 1:

Clear[E1]; eq1 = m0 == E1 + E2 + E3 /. m0 -> K2mass; eq2 = eq1 /. E3 -> (Sqrt[(p3^2) + (m3^2)]) /. m3 -> Kmass; 
    eq3 = eq2 /. p3 -> (Sqrt[(p1^2) + (p2^2) + 2 p1 p2 costh]); 
    eq4 = eq3 /. p1 -> (Sqrt[(E1^2) - (m1^2)]) /. m1 -> Pimass; 
    eq5 = 
      Simplify[eq4 /. p2 -> (Sqrt[(E2^2) - (m2^2)])] /. m2 -> Pimass; 
    solve2 = Solve[eq5, E2];

plot12 = 
      Plot[Evaluate[
          E2 /. solve2[[1]] /. 
            costh -> 1], {E1, Pimass, 
          (((K2mass^2) + (Pimass^2) - (((Kmass + Pimass))^2))/(2.  K2mass
                )) -  .0001}, PlotRange -> {{100, 600}, {100, 600}}]

plot22 = 
      Plot[Evaluate[
          E2 /. solve2[[2]] /. 
            costh -> 1], {E1, Pimass, 
          (((K2mass^2) + (Pimass^2) - (((Kmass + Pimass))^2))/(2.  K2mass
                )) -  .0001}, PlotRange -> {{100, 600}, {100, 600}}]

Show[plot2, plot12, plot22, AspectRatio -> 1]