# Max angle between decay products. A Monte Carlo simulation 
# S.Chekanov.

from java.awt import Color
from jhplot  import *
from java.util import Random
from hephysics.particle import *
from math import *

P=400
# define initial  Lorentz particle
Top=LParticle("Top quark",170)
Top.setV3( 0,0,P )
print Top.toString()

# define Lorentz particles 
W = LParticle("W",80)
b = LParticle("b",5)
q1 = LParticle("q1",0.004)
q2 = LParticle("q2",0.004)

r = Random()
PI=3.1415926
PI2=2*PI
deg= 180 / PI

h1 = H1D("maxAngle",100, 0,180)
for i in range(10000):
     if (i%1000==0): print "event=",i

     Top.twoBodyDecay(W, b,1);
     theta = acos(2.0*r.nextDouble()-1.0 )
     phi   = PI2 *r.nextDouble()
     Top.setThetaPhiP(theta,phi,P)
     W.boost( Top )
     b.boost( Top )
     # W decay 
     W.twoBodyDecay(q1,q2,1)
     # the boost
     q1.boost( W )
     q2.boost( W )
      
     # print q1.toString()
     p=P0D("max angle")     
     p.add(W.angle(q1)*deg)
     p.add(W.angle(q2)*deg)
     pmax=p.getMax()
     h1.fill(pmax)

c1 = HPlot("Canvas",600,400,1, 1)
c1.setGTitle("Max angle t → W b → b q #bar{q}")
c1.visible(1)
c1.setRange(0,180,0,1000)
c1.setNameX("angle_{max} [deg]")
c1.setNameY("Events")

h1.setFill(1)
h1.setFillColor(Color.blue)
c1.draw(h1)

# export to EPS figure 
c1.export("angle"+".eps")