#!/usr/bin/env python

from math import *
import aa
from ana import search
from ROOT import Model, FlatBackground, GaussianPointSource, TFile, TTree, TH1D, TCanvas, AddressOf

# pythonze Model::get_component to work with polymorphism

#wrapped_method = Model.get_component
#def newmethod( name ) :
#    c = wrapped_method(name)
  

bg_only_model = Model()

bg = FlatBackground("background")
bg.norm.value = 1000
bg_only_model.add_component( bg );



sigbg_model = bg_only_model.clone()

g = GaussianPointSource("gauspoint")
g.dec.value   = 0
g.ra.value    = 1
g.norm.value  = 10
g.sigma.value = 1.0 / 180 * pi;
g.sigma.fit   = False

# fix the position
g.ra.fit = g.dec.fit = False;

sigbg_model.add_component( g )

hsig = TH1D("hsig","likelihood ratio for signal experiments",    100, -100,100 )
hbg  = TH1D("hsig","likelihood ratio for background experiments",100, -100,100 )

fout = TFile("outfile.root",'RECREATE')
T    = TTree( 'T', 'pseudoexperiments')
T.Branch('background' , 'FlatBackground'     , AddressOf( bg) )
T.Branch('pointsrc'   , 'GaussianPointSource', AddressOf( g ) )



Ndatasets = 50 # do this many datasets for each Nsig

hists = {}

c = TCanvas('c','test statistic',500,500)




for Nsig in range ( 0, 10 ) :

    H = hists[Nsig] = TH1D("h_likrat"+str(Nsig), "test statistic for Nsig ="+str(Nsig), 100,-1,100 )

    for i in range ( Ndatasets ) :

        # set the true value of norm to the desired number of signal events
        sigbg_model.get_component("gauspoint").norm.value = Nsig

        # generate dataset with exactly Nsig signal events (do Poisson semearing later)
        dataset = sigbg_model.generate_dataset( False )

        # fit the dataset with H0 and H1
        bg_only_model.fit( dataset )
        sigbg_model.fit(   dataset )
        print sigbg_model.info()

        likrat = sigbg_model.fitted_likelihood - bg_only_model.fitted_likelihood
        print ('likelihood ratio', likrat )
        H.Fill( likrat )

        T.Fill()

    if Nsig == 0 : hists[Nsig].SetFillColor(5)
    H.Draw("same" if Nsig>0 else "")
    H.SetLineColor( Nsig )
    c.Update()

fout.Write()