In This Package:

DrawGe68EnergyDiff.py

Go to the documentation of this file.

#!/usr/bin/env python
#
# Draws the net ge68 spectrum and fits the peaks
#
#  Usage:
#   nuwa.py DrawEnergyDiff

# Load DybPython
from DybPython.DybPythonAlg import DybPythonAlg
from GaudiPython import SUCCESS, FAILURE
from GaudiPython import gbl
from decimal import *
from math import sqrt
import array

# Make shortcuts to any ROOT classes you want to use
TF1 = gbl.TF1
TCanvas = gbl.TCanvas
gbl.gStyle.SetPalette(1)
gbl.gStyle.SetOptFit(1)

# Make your algorithm
class DrawEnergyFigsAlg(DybPythonAlg):
    "Algorithm to draw figures from the energy statistics file"
    def __init__(self,name):
        DybPythonAlg.__init__(self,name)
        return

    def initialize(self):
        status = DybPythonAlg.initialize(self)
        if status.isFailure(): return status
        self.info("initializing")

        return SUCCESS

    def execute(self):
        self.info("executing")
        
        return SUCCESS
        
    def finalize(self):
        self.info("finalizing")

        # Make energy histograms
        canvas = TCanvas()

        self.stats["/file0/energy/nobgEnergy"].Draw()
        maxbin = self.stats["/file0/energy/nobgEnergy"].GetMaximumBin()
        max = self.stats["/file0/energy/nobgEnergy"].GetMaximumBin() * 0.006
        interval=0.02
        g = 0
        oldfitChi = 1000
        for i in range(0,15):
          self.stats["/file0/energy/nobgEnergy"].Fit("gaus","","",max-interval,max+2*interval)
          print "interval1=", max-interval
          print "interval2=", max+interval
          fitResults = self.stats["/file0/energy/nobgEnergy"].GetFunction("gaus")
          fitChi = fitResults.GetChisquare()
          ndf = interval*500-3
          fitChi = fitChi / ndf
          print "ndf=", ndf
          print "fitchi=", fitChi
          if oldfitChi > fitChi:
            g = i
          interval=interval+0.02
          oldfitChi = fitChi

        interval=0.02+g*0.02
        print "interval=", interval
        peak1 = TF1("m1","gaus",0,0.75)
        peak2 = TF1("m2","gaus",0.75,1.5)
        total = TF1("mtotal", "gaus(0) + gaus(3)", 0, 1.5)
        #self.stats["/file0/energy/nobgEnergy"].Fit("gaus","","",max-interval,max+2*interval+0.1)
        self.stats["/file0/energy/nobgEnergy"].Fit(peak1,"R")
        self.stats["/file0/energy/nobgEnergy"].Fit(peak2,"R+")
        peak1_0 = peak1.GetParameter(0)
        peak1_1 = peak1.GetParameter(1)
        peak1_2 = peak1.GetParameter(2)


        total = TF1("fit", "[0]*TMath::Gaus(x,[1],[2]) + [3]*TMath::Gaus(x,[1]/2,[2]/sqrt(2))", 0, 1.5)
        total.SetParameter(0,peak1_0)
        total.SetParameter(1,peak1_1)
        total.SetParameter(2,peak1_2)
        total.SetParameter(3,0)

        self.stats["/file0/energy/nobgEnergy"].Fit(total, "R")
        
        # Uncomment to save 
        #f = open('FitData', 'a')
        #fitResults = self.stats["/file0/energy/nobgEnergy"].GetFunction("gaus")
        #for i in range(0,3):
        #  f.write(repr(fitResults.GetParameter(i)))
        #  f.write('\n')
        #f.write('\n')
        #f.close()
        canvas.SaveAs("nobgEnergy.png")
        
        status = DybPythonAlg.finalize(self)
        return status


#####  Job Configuration for nuwa.py ########################################

def configure():
    from StatisticsSvc.StatisticsSvcConf import StatisticsSvc
    statsSvc = StatisticsSvc()
    statsSvc.Input ={"file0":"10000energyStatsDiff_10.root"}
    return

def run(app):
    '''
    Configure and add an algorithm to job
    '''
    app.ExtSvc += ["StatisticsSvc"]
    energyFigsAlg = DrawEnergyFigsAlg("MyEnergyFigs")
    app.addAlgorithm(energyFigsAlg)
    pass

---