<commit_before>import os
import ROOT
import numpy as np
import pyfits
import root_style
import pyLikelihood
from config import get_config
import utils


class Params:
    """Collection of Plotting parameters like Energy bounds,
    colors, file name, etc...."""
    def __init__(self, srcname, Emin=100, Emax=3e5,
                 PlotName="LAT_SED", LineColor=2,
                 PointColor = 1, N = 500):
        self.Emin = Emin #Energy bounds
        self.Emax = Emax
        self.N = N #Number of points for the TGraph
        self.srcname = srcname # Source of interest
        self.PlotName = PlotName #file name
        #color options
        self.LineColor = LineColor
        self.PointColor = PointColor

class Result:
    """Helper class to get the results from a (Un)BinnedAnalysis object
    and compute the SED and errors"""
    def __init__(self, Fit, pars):
        self.Fit = Fit
        self.Model = Fit[pars.srcname].funcs['Spectrum'].genericName()
        self.ptsrc = pyLikelihood.PointSource_cast(Fit[pars.srcname].src)
        self.covar = np.array(utils.GetCovar(pars.srcname, self.Fit, False))
        self.srcpars = pyLikelihood.StringVector()
        Fit[pars.srcname].src.spectrum().getFreeParamNames(self.srcpars)

    def _DumpSED(self,par):
        """Save the energy, E2.dN/dE, and corresponding  error in an ascii file
        The count and residuals plot vs E is also made"""
        E, SED = self.MakeSED(par)
        Err = self.MakeSEDError(par)
        print
        for i in xrange(par.N):
          if Err[i]/SED[i] == min(Err/SED):
            print "Decorrelation energy : %4.2e MeV"%E[i]
            print "Diffential flux  at the Decorrelation energy : %2.2e +/-  %2.2e ph/cm2/s/MeV" %(SED[i]/E[i]**2*624152.206,Err[i]/E[i]**2*624152.206)
            print "SED value at the Decorrelation energy : %2.2e +/-  %2.2e erg/cm2/s" %(SED[i],Err[i])
        print

        try:
            self.CountsPlot(par)
        except:
            pass
        # Save all in ascii file
        # log(E)  log (E**2*dN/dE)   log(E**2*dN/dE_err)  is_dot (0,1) is_upper (0,1)
        save_file = open(par.PlotName + '.dat', 'w')
        save_file.write("log(E)  log (E**2*dN/dE)   Error on log(E**2*dN/dE)   \n")
        for i in xrange(par.N):
            save_file.write("%12.4e  %12.4e  %12.4e \n" % (E[i], SED[i], Err[i]))
        save_file.close()

    def MakeFlux(self, params):
        """Compute differential Flux distribution and
        corresponding energy and return a numpy array"""
        E = np.logspace(np.log10(params.Emin), np.log10(params.Emax), params.N)
        Flux = np.zeros(params.N)
        for i in xrange(params.N):
            Flux[i] = self.dNde(E[i])
        return E, Flux

    def MakeSED(self, pars):
        """Compute Spectral energy distribution and corresponding energy
        and return a numpy array"""
        E = np.logspace(np.log10(pars.Emin), np.log10(pars.Emax), pars.N)
        nuFnu = np.zeros(pars.N)
        for i in xrange(pars.N):
            nuFnu[i] = 1.602e-6 * E[i] ** 2 * self.dNde(E[i]) #Mev to Ergs
        return E, nuFnu

    def MakeSEDError(self, pars):
        """@todo: document me"""
        estep = np.log(pars.Emax / pars.Emin) / (pars.N - 1)
        energies = pars.Emin * np.exp(estep * np.arange(np.float(pars.N)))
        err = np.zeros(pars.N)
        j = 0
        for ene in energies:
            arg = pyLikelihood.dArg(ene)
            partials = np.zeros(len(self.srcpars))
            for i in xrange(len(self.srcpars)):
                x = self.srcpars[i]
                partials[i] = self.ptsrc.spectrum().derivByParam(arg, x)
            err[j] = np.sqrt(np.dot(partials, np.dot(self.covar, partials)))
            j += 1

        return 1.602e-6 * energies ** 2 * err #Mev to Ergs

    def dNde(self, energy):
        arg = pyLikelihood.dArg(energy)
        return self.ptsrc.spectrum()(arg)

    def CountsPlot(self, Parameter):
        """@todo: document me"""
        ROOT.gROOT.SetBatch(ROOT.kTRUE)
        imName = "tmp.fits"
        self.Fit.writeCountsSpectra(imName)
        image = pyfits.open(imName)

        #loop on the source names to find the good one
        j = 0
        for ID in image[1].data.names:
            if ID == Parameter.srcname:
                indice = j
            j += 1

        emax = image[3].data.field(0)
        emin = image[3].data.field(1)

        E = np.array((emax + emin) / 2.)
        err_E = np.array((-emax + emin) / 2.)

        src = np.array(image[1].data.field(indice))
        Nbin = len(src)

        obs = np.array(image[1].data.field(0))
        obs_err = np.array(np.sqrt(image[1].data.field(0)))
        total = 0
        for i in xrange(len(image[1].data.names) - 1):
            total = total + image[1].data.field(i + 1)

        other = np.array(total - image[1].data.field(indice))
        total = np.array(total)
        residual = np.zeros(Nbin)
        Dres = np.zeros(Nbin)

        for i in xrange(Nbin):
            try:
                residual[i] = (obs[i] - total[i]) / total[i]
                Dres[i] = (obs_err[i] / total[i])
            except:
                residual[i] = 0.
                Dres[i] = 0.

        cplot = ROOT.TCanvas("Count Plot")
        cplot.SetLogy()
        cplot.SetLogx()
        ghcount = ROOT.TH2F("ghcount", "", 80, 91, 499e3, 100, 0.1, max(obs) * 2)
        ghcount.SetStats(000)
        ghcount.SetXTitle("E (MeV) ")
        ghcount.SetYTitle("Counts / bin")
        ghcount.Draw()
        tgrobs = ROOT.TGraphErrors(Nbin, E, obs, err_E, obs_err)
        tgrobs.SetLineColor(1)
        tgrobs.SetMarkerColor(1)
        tgrobs.SetMarkerStyle(1)
        tgrobs.Draw("P")

        tgrother = ROOT.TGraph(Nbin, E, other)
        tgrother.SetLineWidth(2)
        tgrother.SetLineStyle(2)
        tgrother.Draw("L")

        tgr = ROOT.TGraph(Nbin, E, src)
        tgr.SetLineColor(1)
        tgr.SetLineWidth(2)
        tgr.Draw("L")

        tgrsum = ROOT.TGraph(Nbin, E, total)
        tgrsum.SetLineStyle(3)
        tgrsum.SetLineWidth(2)
        tgrsum.Draw("L")

        cres = ROOT.TCanvas("Residuals plot")
        cres.SetLogx()
        ymin = min(residual) - max(Dres)
        ymax = max(residual) + max(Dres)
        ghres = ROOT.TH2F("ghres", "", 80, 91, 499e3, 100, ymin, ymax)
        ghres.SetStats(000)
        ghres.SetXTitle("E (MeV) ")
        ghres.SetYTitle("(counts -model)/model")
        ghres.Draw()
        tgres = ROOT.TGraphErrors(Nbin, E, residual, err_E, Dres)
        tgres.Draw("P")

        zero = np.zeros(2)
        Ezero = np.array([0, 1e10])
        tg0 = ROOT.TGraph(2, Ezero, zero)
        tg0.SetLineStyle(2)
        tg0.Draw("L")

        # Save the plots in different formats
        filebase = Parameter.PlotName
        cplot.Print(filebase + "_CountsPlot.eps")
        cplot.Print(filebase + "_CountsPlot.C")
        cploy.Print(filebase + "_CountsPlot.png")
        cres.Print(filebase + "_ResPlot.eps")
        cres.Print(filebase + "_ResPlot.C")
        cres.Print(filebase + "_ResPlot.png")
        os.system("rm " + imName)
        image.close()

def PlotTS(Time, TimeErr, TS):
    """Scatter plot TS(Time)"""
    Time = np.asarray(Time)
    TimeErr = np.asarray(TimeErr)
    TS = np.asarray(TS)
    zero = np.zeros_like(TS)
    gh = ROOT.TH2F("ghts", "", 80, min(Time) - max(TimeErr) * 3,
                   max(Time) + max(TimeErr) * 3, 100, 0, max(TS) * 1.2)
    gh.SetStats(000)<commit_msg>fix res and count plot<commit_after>import os
import ROOT
import numpy as np
import pyfits
import root_style
import pyLikelihood
from config import get_config
import utils
import array as ar

class Params:
    """Collection of Plotting parameters like Energy bounds,
    colors, file name, etc...."""
    def __init__(self, srcname, Emin=100, Emax=3e5,
                 PlotName="LAT_SED", LineColor=2,
                 PointColor = 1, N = 500):
        self.Emin = Emin #Energy bounds
        self.Emax = Emax
        self.N = N #Number of points for the TGraph
        self.srcname = srcname # Source of interest
        self.PlotName = PlotName #file name
        #color options
        self.LineColor = LineColor
        self.PointColor = PointColor

class Result:
    """Helper class to get the results from a (Un)BinnedAnalysis object
    and compute the SED and errors"""
    def __init__(self, Fit, pars):
        self.Fit = Fit
        self.Model = Fit[pars.srcname].funcs['Spectrum'].genericName()
        self.ptsrc = pyLikelihood.PointSource_cast(Fit[pars.srcname].src)
        self.covar = np.array(utils.GetCovar(pars.srcname, self.Fit, False))
        self.srcpars = pyLikelihood.StringVector()
        Fit[pars.srcname].src.spectrum().getFreeParamNames(self.srcpars)

    def _DumpSED(self,par):
        """Save the energy, E2.dN/dE, and corresponding  error in an ascii file
        The count and residuals plot vs E is also made"""
        E, SED = self.MakeSED(par)
        Err = self.MakeSEDError(par)
        print
        for i in xrange(par.N):
          if Err[i]/SED[i] == min(Err/SED):
            print "Decorrelation energy : %4.2e MeV"%E[i]
            print "Diffential flux  at the Decorrelation energy : %2.2e +/-  %2.2e ph/cm2/s/MeV" %(SED[i]/E[i]**2*624152.206,Err[i]/E[i]**2*624152.206)
            print "SED value at the Decorrelation energy : %2.2e +/-  %2.2e erg/cm2/s" %(SED[i],Err[i])
        print

        try:
            self.CountsPlot(par)
        except:
            pass
        # Save all in ascii file
        # log(E)  log (E**2*dN/dE)   log(E**2*dN/dE_err)  is_dot (0,1) is_upper (0,1)
        save_file = open(par.PlotName + '.dat', 'w')
        save_file.write("log(E)  log (E**2*dN/dE)   Error on log(E**2*dN/dE)   \n")
        for i in xrange(par.N):
            save_file.write("%12.4e  %12.4e  %12.4e \n" % (E[i], SED[i], Err[i]))
        save_file.close()

    def MakeFlux(self, params):
        """Compute differential Flux distribution and
        corresponding energy and return a numpy array"""
        E = np.logspace(np.log10(params.Emin), np.log10(params.Emax), params.N)
        Flux = np.zeros(params.N)
        for i in xrange(params.N):
            Flux[i] = self.dNde(E[i])
        return E, Flux

    def MakeSED(self, pars):
        """Compute Spectral energy distribution and corresponding energy
        and return a numpy array"""
        E = np.logspace(np.log10(pars.Emin), np.log10(pars.Emax), pars.N)
        nuFnu = np.zeros(pars.N)
        for i in xrange(pars.N):
            nuFnu[i] = 1.602e-6 * E[i] ** 2 * self.dNde(E[i]) #Mev to Ergs
        return E, nuFnu

    def MakeSEDError(self, pars):
        """@todo: document me"""
        estep = np.log(pars.Emax / pars.Emin) / (pars.N - 1)
        energies = pars.Emin * np.exp(estep * np.arange(np.float(pars.N)))
        err = np.zeros(pars.N)
        j = 0
        for ene in energies:
            arg = pyLikelihood.dArg(ene)
            partials = np.zeros(len(self.srcpars))
            for i in xrange(len(self.srcpars)):
                x = self.srcpars[i]
                partials[i] = self.ptsrc.spectrum().derivByParam(arg, x)
            err[j] = np.sqrt(np.dot(partials, np.dot(self.covar, partials)))
            j += 1

        return 1.602e-6 * energies ** 2 * err #Mev to Ergs

    def dNde(self, energy):
        arg = pyLikelihood.dArg(energy)
        return self.ptsrc.spectrum()(arg)

    def CountsPlot(self, Parameter):
        """@todo: document me"""
        ROOT.gROOT.SetBatch(ROOT.kTRUE)
        imName = "tmp.fits"
        self.Fit.writeCountsSpectra(imName)
        image = pyfits.open(imName)

        #loop on the source names to find the good one
        j = 0
        for ID in image[1].data.names:
            if ID == Parameter.srcname:
                indice = j
            j += 1

        emax = image[3].data.field(0)
        emin = image[3].data.field(1)

        E = np.array((emax + emin) / 2.)
        err_E = np.array((emax - emin) / 2.)

        src = np.array(image[1].data.field(indice))
        Nbin = len(src)

        obs = np.array(image[1].data.field(0))
        obs_err = np.array(np.sqrt(image[1].data.field(0)))
        total = 0
        for i in xrange(len(image[1].data.names) - 1):
            total = total + image[1].data.field(i + 1)

        other = np.array(total - image[1].data.field(indice))
        total = np.array(total)
        residual = np.zeros(Nbin)
        Dres = np.zeros(Nbin)

        cplot = ROOT.TCanvas("Counts_Plot")
        cplot.SetLogy()
        cplot.SetLogx()
        ghcount = ROOT.TH2F("ghcount", "", 80, min(E)*0.3,max(E)*2, 100, 0.1, max(obs) * 2)
        ghcount.SetStats(000)
        ghcount.SetXTitle("E (MeV) ")
        ghcount.SetYTitle("Counts / bin")
        ghcount.Draw()

        tgrobs = ROOT.TGraphErrors(Nbin, ar.array('f',E), ar.array('f',obs), ar.array('f',err_E), ar.array('f',obs_err))
        tgrobs.SetLineColor(2)
        tgrobs.SetMarkerColor(2)
        tgrobs.SetMarkerStyle(20)
        tgrobs.Draw("pz")

        tgrother = ROOT.TGraph(Nbin, ar.array('f',E), ar.array('f',other))
        tgrother.SetLineWidth(2)
        tgrother.SetLineStyle(2)
        tgrother.Draw("L")

        tgr = ROOT.TGraph(Nbin, ar.array('f',E), ar.array('f',src) )
        tgr.SetLineColor(1)
        tgr.SetLineWidth(2)
        tgr.Draw("L")

        tgrsum = ROOT.TGraph(Nbin,  ar.array('f',E),  ar.array('f',total))
        tgrsum.SetLineStyle(3)
        tgrsum.SetLineWidth(2)
        tgrsum.Draw("L")

        legarc =ROOT.TLegend(0.6968391,0.6716102,0.8807471,0.845339);
        legarc.AddEntry(tgrobs,"Data","lp");
        legarc.AddEntry(tgr,Parameter.srcname,"l");
        legarc.AddEntry(tgrother,"Other Sources","l");
        legarc.AddEntry(tgrsum,"All Sources","l");
        legarc.SetFillColor(0)
        legarc.Draw()

        cres = ROOT.TCanvas("Residuals_Plot")
        cres.SetLogx()

        for i in xrange(Nbin):
            try:
                residual[i] = (obs[i] - total[i]) / total[i]
                Dres[i] = (obs_err[i] / total[i])
            except:
                residual[i] = 0.
                Dres[i] = 0.
            if residual[i] == -1.:
               residual[i] = 0.

        ymin = min(residual) - max(Dres)
        ymax = max(residual) + max(Dres)
        ghres = ROOT.TH2F("ghres", "", 80, min(E)*0.3,max(E)*2, 100, ymin, ymax)
        ghres.SetStats(000)
        ghres.SetXTitle("E (MeV) ")
        ghres.SetYTitle("(counts -model)/model")
        ghres.Draw()
        tgres = ROOT.TGraphErrors(Nbin, ar.array('f',E), ar.array('f',residual), ar.array('f',err_E), ar.array('f',Dres))
        tgres.SetLineColor(2)
        tgres.SetMarkerColor(2)
        tgres.Draw("p*z")

        zero = np.zeros(2)
        Ezero = np.array([0, 1e10])
        tg0 = ROOT.TGraph(2, ar.array('f',Ezero), ar.array('f',zero))
        tg0.SetLineStyle(2)
        tg0.Draw("L")

        # Save the plots in different formats
        filebase = Parameter.PlotName
        cplot.Print(filebase + "_CountsPlot.eps")
        cplot.Print(filebase + "_CountsPlot.C")
        cplot.Print(filebase + "_CountsPlot.png")
        cres.Print(filebase + "_ResPlot.eps")
        cres.Print(filebase + "_ResPlot.C")
        cres.Print(filebase + "_ResPlot.png")
        os.system("rm " + imName)
        image.close()

def PlotTS(Time, TimeErr, TS):
    """Scatter plot TS(Time)"""
    Time = np.asarray(Time)
    TimeErr = np.asarray(TimeErr)
    TS = np.asarray(TS)
    zero = np.zeros_like(TS)
    gh = ROOT.TH2F("ghts", "", 80, min(Time) - max(TimeErr) * 3,
                   max(Time) + max(TimeErr) * 3, 100, 0, max(TS) * 1.2)
    gh.SetStats(000)
{'fraction_non_alphanumeric': 0.1453941484604574, 'fraction_numerical': 0.023827775648396575, 'mean_word_length': 2.096933728981207, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '6769983', 'n_tokens_mistral': 5983, 'n_tokens_neox': 5617, 'n_words': 1164}